Understanding Kernel Ridge Regression:
Common behaviors from simple functions to density functionals

Machine learning (ML) is a powerful data-driven method for learning patterns in high-dimensional spaces via induction, and has had widespread success in many fields including more recent applications in quantum chemistry and materials science Müller et al. (2001); Kononenko (2001); Sebastiani (2002); Ivanciuc (2007); Bartók et al. (2010); Rupp et al. (2012); Pozun et al. (2012); Hansen et al. (2013); Schütt et al. (2014). Here we are interested in applications of ML to construction of density functionals Snyder et al. (2012, 2013a); Li et al. (2014); Snyder et al. (2013b, 2015), which have focused so far on approximating the kinetic energy (KE) of non-interacting electrons. An accurate, general approximation to this could make orbital-free DFT a practical reality.

However, ML methods have been developed within the areas of statistics and computer science, and have been applied to a huge variety of data, including neuroscience, image and text processing, and robotics Burke (2006). Thus, they are quite general and have not been tailored to account for specific details of the quantum problem. For example, it was found that a standard method, kernel ridge regression, could yield excellent results for the KE functional, while never yielding accurate functional derivatives. The development of methods for bypassing this difficulty has been important for ML in general Snyder et al. (2013b).

ML provides a whole suite of tools for analyzing data, fitting highly non-linear functions, and dimensionality reduction Hastie et al. (2009). More importantly in this context, ML provides a completely different way of thinking about electronic structure. The traditional ab-initio approach Levine (2009) to electronic structure involves deriving carefully constructed approximations to solving the Schrödinger equation, based on physical intuition, exact conditions and asymptotic behaviors Dreizler and Gross (1990). On the other hand, ML learns by example. Given a set of training data, ML algorithms learn via induction to predict new data. ML provides limited interpolation over a specific class of systems for which training data is available.

A system of $N$ interacting electrons with some external potential is characterized by a $3N$ coordinate wavefunction, which becomes computationally demanding for large $N$. In the mid 1960’s, Hohenberg and Kohn proved a one-to-one correspondence between the external potential of a quantum system and its one-electron ground-state density Hohenberg and Kohn (1964), showing that all properties are functionals of the ground-state density alone, which can in principle be found from a single Euler equation for the density. Although these fundamental theorems of DFT proved the existence of a universal functional, essentially all modern calculations use the KS scheme Kohn and Sham (1965), which is much more accurate, because the non-interacting KE is found exactly by using an orbital-scheme Pribram-Jones et al. (2014). This is far faster than traditional approaches for large $N$, but remains a bottleneck. If a sufficiently accurate density functional for the non-interacting electrons could be found, it could increase the size of computationally tractable systems by orders of magnitude.

The Hohenberg-Kohn theorem guarantees that all properties of the system can be determined from the electronic density alone. The basic tenet of ML is that a pattern must exist in the data in order for learning to be possible. Thus, DFT seems an ideal case to apply ML. ML learns the underlying pattern in solutions to the Schrödinger equation, bypassing the need to directly solve it. The HK theorem is a statement concerning the minimal information needed to do this for an arbitrary one-body potential.

Some of us recently used ML to learn the non-interacting KE of fermions in a one-dimensional box subject to smooth external potentials Snyder et al. (2012) and of a one-dimensional model of diatomics where we demonstrated the ability of ML to break multiple bonds self-consistently via an orbital-free density functional theory (DFT) Snyder et al. (2013a). Such KE data is effectively noise-free, since it is generated via deterministic reference calculations, by solving the Schrödinger equation or KS equations numerically exactly. (The limited precision of the calculation might be considered “noise,” as different implementations might yield answers differing on the order of machine precision, but this is negligble.) There is no noise, in the traditional sense, as is typically associated with experimental data. Note that what is considered “noise” depends on what is considered ground truth, i.e., the data to be learned. In particular, if a single reference method is used, its error with respect to a universal functional is not considered noise for the ML model. A perfect ML model should, at best, precisely reproduce the single-reference calculation.

As an example, in Fig. 1 we plot a measure of the error of ML for the KE of up to 4 noninteracting spinless fermions in a box under a potential with 9 parameters (given in detail in Ref. Snyder et al. (2012)), fitted for different numbers of evenly spaced training densities as a function of the hyperparameter $\sigma$ (called the length scale), for fixed $\lambda$ (a hyperparameter called the regularization strength) and several different number of training points $N_{T}$. The scale is logarithmic ¹¹1 We will use $\log$ to denote $\log_{10}$ here and throughout this work, so there are large variations in the fitted error. We will give a more in-depth analysis of the model performance on this data set in a later section after we have formally defined the functions and hyperparameters involved, but for now it is still useful to observe the qualitative behaviors that emerge in the figure. Note that the curves assume roughly the same shape for each $N_{T}$ over the range of $\sigma$, and that they all possess distinct features in different regimes of $\sigma$.

To better understand the behavior with respect to hyperparameters seen in Fig 1, we have chosen in this paper to apply them to the prototypical regression problem, that of fitting a simple function of one coordinate. We also remove all stochastic elements of the procedure, by considering data points on uniform grids, defining errors in the continuum limit, etc. This is shown in Fig. 2, where we plot a measure of the error of ML for a simple function $\cos x$, fitted in the region between 0 and 1, inclusive, for several $N_{T}$ (represented as values on a grid) as a function of $\sigma$. Note the remarkable similarity between the features and characteristics of the curves of this figure and those of Fig. 1 (like Fig. 1 before it, we will give a more in-depth analysis of Fig. 2 later). We explore the behavior of the fitting error as a function of the number of training parameters and the hyperparameters that are used in kernel ridge regression with Gaussian kernel. We find the landscape to be surprisingly rich, and we also find elegant simplicities in various limiting cases. After this, we will be able to characterize the behavior of ML for systems like the one shown in Fig. 1.

Looking at Fig. 2, we see that the best results (lowest error) are always obtained from the middle of the curves, which can become quite flat with enough training data. Thus, any method for determining hyperparameters should usually yield a length scale somewhere in this valley. For very small length scales, all curves converge to the same poor result, regardless of the number of training points. On the other hand, notice also the plateau structure that develops for very large length scales, again with all curves converging to a certain limit. We show for which ranges of hyperparameters these plateaus emerge and how they can be estimated. We also study and explain many of the features of these curves. To show the value of this study, we then apply the same reasoning to the problem that was tackled in Ref. Snyder et al. (2012); Li et al. (2014), which we showcased in Fig. 1. From the machine learning perspective our study may appear unusual as it considers properties in data and problems that are uncommon. Namely there are only a few noise free data points and all are low dimensional. Nevertheless, from the physics point of view the toy data considered reflects very well the essential properties of a highly relevant problem in quantum physics: the machine learning of DFT.

In this work, we will first use ML to fit a very simple function of one variable,

on the interval $x\in[0,1]$. We will focus on exploring the properties of ML for this simple function before proceeding to our DFT cases. We choose a set of $x$-values and corresponding $f(x)$ values as the “training data” for ML to learn from. In ML, the $x$-values $\{x_{j}\}$ for $j=1,\dots,N_{T}$ are known as features, and corresponding $f$-values, $\{f_{j}\}$, are known as labels. Here $N_{T}$ is the number of training samples. Usually, ML is applied with considerable random elements, such as in the choice of data points and selection of test data. In our case, we choose evenly spaced training points on the interval $[0,1]$: $x_{j}=(j-1)/(N_{T}-1)$ for $j=1,\dots,N_{T}$.

Using this dataset, we apply kernel ridge regression (KRR), which is a non-linear form of regression with regularization to prevent overfitting Hastie et al. (2009), to fit $f(x)$. The general form of KRR is

where $\alpha_{j}$ are the weights, $k$ is the kernel (which is a measure of similarity between features), and the hyperparameter $\lambda$ controls the strength of the regularization and is linked to the noise level of the learning problem. We use the Gaussian kernel

a standard choice in ML that works well for a variety of problems. The hyperparameter $\sigma$ is the length scale of the Gaussian, which controls the degree of correlation between training points.

The weights $\alpha_{j}$ are obtained through the minimization of the cost function

where

The exact solution is given by

where $\mathbfit{I}$ is the $N_{T}\times N_{T}$ identity matrix, $\mathbfit{K}$ is the kernel matrix with elements $K_{ij}$ = $k(x_{i},x_{j})$, and $\mathbfit{f}=(f_{1},\dots,f_{N_{T}})^{T}$.

The two parameters $\lambda$ and $\sigma$ not determined by Eq. (6) are called hyperparameters and must be determined from the data (see section III.3). $\sigma$ can be viewed as the characteristic length scale of the problem being learned (the scale on which changes of $f$ take place), as discernible from the data (and thus dependent on, e.g., the number of training samples). $\lambda$ controls the leeway the model has to fit the training points. For small $\lambda$, the model has to fit the training points exactly, whereas for larger $\lambda$ some deviation is allowed. Larger values of $\lambda$ therefore cause the model to be smoother and vary less, i.e., less prone to overfitting. This can be directly seen in Gaussian process regression Rasmussen and Williams (2006), a related Bayesian ML model with predictions identical to those of KRR. There, $\lambda$ formally is the variance of the assumed additive Gaussian noise in values of $f$.

KRR is a method of interpolation. Here, we are mainly concerned with the error of the machine learning approximation (MLA) to $f(x)$ in the interpolation region, which in this case is the interval $x\in[0,1]$. As a measure of this error, we define

In the case of the Gaussian kernel, we can expand this and derive the integrals that appear analytically. To simplify the analytical process, we define

as the benchmark error when $f^{\rm ML}(x)\equiv 0$. For the cosine function in Eq. (1),

Our goal is to characterize the dependence of the performance of the model on the size of the training data set $(N_{T})$ and the hyperparameters of the model $(\sigma,\,\lambda)$. For this simple model, we discuss different regions of qualitative behavior and derive the dependence of $\Delta f$ for various limiting values of these hyperparameters; we do all of this in the next few sections. In Section IV, we discuss how these results can be qualitatively generalized for the problem of using ML to learn the KE functional for non-interacting fermions in the box for a limited class of potentials.

We begin by analyzing the structure of $\Delta f$ as a function of $\sigma$ for fixed $\lambda$ and $N_{T}$. Fig. 2 shows $\Delta f$ as a function of $\sigma$ for various $N_{T}$ while fixing $\lambda=10^{-6}$. The curves have roughly the same “valley” shape for all $N_{T}$. The bottom of the valley is an order of magnitude deeper than the walls and may have multiple local minima. These valleys are nearly identical in shape for sufficiently large $N_{T}$, which indicates that this particular feature arises in a systematic manner as $N_{T}$ increases. Moreover, this central valley opens up to the left (i.e., smaller $\sigma$) as $N_{T}$ increases— as the training samples become more densely packed, narrower Gaussians are better able to interpolate the function smoothly. Thus, with more training samples, a wider range of $\sigma$ values will yield an accurate model.

In addition, a “cusp” will begin to appear in the region to the left of the valley, and its general shape remains the same for increasing $N_{T}$. This is another recurring feature that appears to develop systematically like the valley. For a fixed $N_{T}$, and starting from the far left, the $\Delta f$ curve begins to decrease monotonically to the right, i.e., as $\sigma$ increases. The cusps mark the first break in this monotonic behavior, as $\Delta f$ increases briefly after reaching this local minimum before resuming its monotonic decrease for increasing $\sigma$ (until this monotonicity is interrupted again in the valley region). The cusps shift to the left as $N_{T}$ increases, following the trend of the valleys. This indicates that they are a fundamental feature of the $\Delta f$ curves and that their appearances coincide with a particular behavior of the model as it approaches certain $\sigma$ values. Note that $\Delta f$ decreases nearly monotonically as $N_{T}$ increases for all $\sigma$. This is as expected, since each additional training sample adds another weighted Gaussian, which should improve the fit.

Fig. 3 again shows $\Delta f$ as a function of $\sigma$, but for various $\lambda$ with $N_{T}$ fixed at 33. As $\lambda$ decreases, $\Delta f$ again decreases nearly monotonically and the central region opens up to the right (i.e., larger $\sigma$). Note that the curves for each $\lambda$ coincide up to a certain $\sigma$ before they split off from the rest, with the lower $\lambda$-valued curves breaking off further along to the right than those with larger $\lambda$. This shows a well-known phenomenon, namely that regularization strength $\lambda$ and kernel length scale $\sigma$ both give rise to regularization and smoothing Smola et al. (1998). Additionally, we observe the emergence of “plateau”-like structures on the right. These will be explored in detail in Section III.1.2.

A canonical quantum system used frequently to explore basic quantum principles and as a proving ground for approximate quantum methods is the particle in a box. In this case, we confine one fermion to a 1d box with hard walls at $x=0,1$, with the addition of the external potential $v(x)$ in the interval $x\in[0,1]$. The equation that governs the quantum mechanics is the familiar one-body Schrödinger equation

A solution of this equation gives the orbitals $\phi_{j}(x)$ and energies $\epsilon_{j}$. For one fermion, only the lowest energy level is occupied. The total energy is $E=\epsilon_{1}$, the potential energy is

(where $n(x)=|\phi(x)|^{2}$ is the electron density), and the KE is $T=E-V$. In the case of one particle, the KE can be expressed exactly in terms of the electron density, known as the von Weizsäcker functional Weizsäcker (1935)

where $n^{\prime}(x)=dn/dx$. Our goal here in this section is not to demonstrate the efficacy of ML approximations for the KE in DFT (which is the subject of other works Snyder et al. (2012, 2013a)), but rather to study the properties of the ML approximations with respect to those applications.

We choose a simple potential inside the box,

to model a well of depth $D$, which has also been used in the study of semiclassical methods Cangi et al. (2011). To generate reference data for ML to learn from, we solve Eq. (41) numerically by discretizing space onto a uniform grid, $x_{j}=(j-1)/(N_{G}-1)$, for $j=1,\dots,N_{G}$, where $N_{G}$ is the number of grid points. Numerov’s method is used to solve for the lowest energy orbital and its corresponding eigenvalue. We compute $T$ and $n(x)$, which is represented by its values on the grid. For a desired number of training samples $N_{T}$, we vary $D$ uniformly over the range $[0,100]$, inclusive, generating $N_{T}$ pairs of electron densities and exact KEs. Additionally, a test set with $500$ pairs of electron densities and exact KEs is generated.

As in the previous sections, we use KRR to learn the KE of this model system. The formulation is identical to that of Ref. Snyder et al. (2012):

where $k$ is the Gaussian kernel

and

where $\Delta x=1/(N_{G}-1)$ is the grid spacing. The weights are again given by Eq. (6), found by minimizing the cost function in Eq. (4).

In analogy to Eq. (7), we measure the error of the model as the total squared error integrated over the interpolation region

where $n_{D}$ is the exact density for the potential with well depth $D$, and $T[n_{D}]$ is the exact corresponding KE. We approximate the integral by Simpson’s rule evaluated on $D$ sampled over the test set (i.e., 500 values spaced uniformly over the interpolation region). This sampling is sufficiently dense over the interval to give an accurate approximation to $\Delta T$.

In Fig. 17, we plot $\Delta T$ as a function of the length scale of the Gaussian kernel, $\sigma$, for various training set sizes $N_{T}$. Clearly, the trends are very similar to Fig. 2: the transition $\sigma_{s}$ between the regions I and II becomes smaller as $N_{T}$ increases, the valley in region II widens, and region III on the right remains largely unchanged. The dependence of $\sigma_{s}$ on $N_{T}$ appears to follow the same power law $\sigma_{s}\propto N_{T}^{p}$, but the value of $p$ is different in this case. A rough estimate yields $p\approx-0.8$, which is similar to $p=-1$ for the simple cosine function explored in the previous sections, but the details will depend on the specifics of the data.

Similarly, Fig. 18 shows the same plot but with $N_{T}$ fixed and $\lambda$ varied. Again, the same features are present as in Fig. 3, i.e., three regions with different qualitative behaviors. In region I, $\Delta T$ has the same decay shape as the kernel functions (Gaussians) begin to overlap significantly, making it possible for the regression to function properly and fit the data. For large values of the regularization strength $\lambda$, the model over-regularizes, yielding high errors for any value of $\sigma$. As $\lambda$ decreases, the weights are given more flexibility to conform to the shape of KE functional. Using the same definition for the estimation of $\sigma_{s}$ in Eq. (10), the median nearest neighbor distance over this training set gives $\sigma_{s}=0.019$. We then have $\log\sigma_{s}=-1.72$, which matches the boundary between regions I and II in Fig. 18. In region III, the same plateau features emerge for small $\lambda$. Again, these plateaus occur when polynomial forms of the regression model become dominant for certain combinations of the parameters $\sigma$, $\lambda$, and $N_{T}$.

From Eq. (12) and Eq. (13), we showed that the model error will tend to the benchmark error while $\sigma\to 0$ or $\lambda\to\infty$. Similarly to Eq. (8), we can also define the benchmark error when $T^{\rm ML}[n]\equiv 0$ for this data set as

Evaluating the above integral numerically on the test set, we obtain $\log\Delta T_{0}=3.7$. This matches the error when $\sigma\to 0$ in Fig. 17 and Fig. 18.

We define the $\sigma$ that gives the global minimum of $\Delta T$ as $\tilde{\sigma}$; similarly to Fig. 12, we plot the optimal model error $\Delta T(\tilde{\sigma})$ as a function of $\lambda$ in Fig. 19.

For large $\lambda$, we overregularize the model; the model error tends to the benchmark error in Eq. (49). For moderate $\lambda$, we observe the same linear proportionality $\Delta T(\tilde{\sigma})\propto\lambda$ as in Fig. 12.

In this toy system, the prediction of the KE from KRR models shares properties similar to those that we explored in learning the 1d cosine function. Now we will consider up to 4 noninteracting spinless fermions under a potential with 9 parameters as reported in Ref. Snyder et al. (2012).

These densities are represented on $N_{G}=500$ evenly spaced grid points in $0\leq x\leq 1$. Here a model is built using $N_{T}/4$ pairs of electron densities and exact KEs for each particle number $N=1,2,3,4$, respectively. Thus, the size of the training set is $N_{T}$. 1000 pairs of electron densities and exact KEs are generated for each $N$, so the size of the test set is $S=4000$. Since there are 9 parameters in the potential, we cannot define the error as an integral, so we use summation instead. Thus, the error on the test set is defined as the mean square error (MSE) on the test densities

Fig. 1 shows the error of the model as a function of $\sigma$ with various $N_{T}$ for fixed $\lambda=10^{-10}$. Although this system is more complicated than the previous two systems discussed in this paper, the qualitative behaviors in Fig. 1 are similar to Fig. 2 and Fig. 17. Table I in Ref. Snyder et al. (2012) only shows the model error with optimized hyperparameters for $N_{T}=400$. In Fig. 20, model errors varying with a wide range of $\sigma$ values are shown for 4 different values of $\lambda$ ²²2For large $\sigma$, if $\lambda$ is small, there will be numerical instability in the computation of $(\mathbfit{K}+\lambda\mathbfit{I})^{-1}$. Thus, the $\lambda=3.2\times 10^{-14}$ curve is not plotted for $\log\sigma>5.5$.. The qualitative properties in Fig. 20 are similar to Fig. 3 and Fig. 18. For example, the existence of three regions with distinctly different behavior for the model error can be ascertained just like before. In region I, error curves with different $\lambda$ will all tend to the same benchmark error limit when $\sigma\to 0$. The median nearest neighbor distance over this training set gives $\sigma_{s}=0.022$. In Fig. 20, the boundary between region I and region II is well estimated by $\log\sigma_{s}=-1.66$. In region III, the familiar plateau features emerge. In region II, where $\sigma$ is such that the model is optimal or close to it, we find that the model with hyperparameters $\sigma=1.86,\,\lambda=3.2\times 10^{-14}$ performs the best. The MSE for this model is $1.43\times 10^{-7}$ Hartree. Another common measure of error is the mean absolute error (MAE), which is also used in Ref. Snyder et al. (2012). The MAE of this model is $1.99\times 10^{-4}\,\text{Hartree}=0.12\,\text{kcal}/\text{mol}$. This result is consistent with the model performance reported in Ref. Snyder et al. (2012)³³3 In Ref. Snyder et al. (2012), the densities were treated as vectors. Here, we treated the densities as functions, so the length scale mentioned here is related to the length scale in Ref. Snyder et al. (2012) by a factor of $\sqrt{\Delta x}=0.045$..

In this work, we have analyzed the properties of KRR models with a Gaussian kernel applied to fitting a simple 1d function. In particular, we have explored regimes of distinct qualitative behavior and derived the asymptotic behavior in certain limits. Finally, we generalized our findings to the problem of learning the KE functional of a single particle confined to a box and subject to a well potential with variable depth. Considering the vast difference in nature of the two problems compared in this work, a 1d cosine function and the KE as a functional of the electron density (a very high-dimensional object), the similarities of the measures of error $\Delta f$ and $\Delta T$ between each other are remarkable. This analysis demonstrates that much of the behavior of the model can be rationalized by and distilled down to the properties of the kernel. Our goal in this work was to deepen our understanding of how the performance of KRR depends on the parameters qualitatively, in particular in the case that is relevant for MLA in DFT, namely the one of noise-free data and high-dimensional inputs, and how one may determine a-priori which regimes the model lies in. From the ML perspective the scenario analyzed in this work was an unusual one: small data, virtually no noise, low dimensions and high complexity. The effects found are not only interesting from the physics perspective, but are also illuminating from a learning theory point of view. However, in ML practice the extremes that contain plateaus or the “comb” region will not be observable, as the practical data with its noisy manifold structure will confine learning in the favorable region II. Future work will focus on theory and practice in order to improve learning techniques for low noise problems.