{'fraction_non_alphanumeric': 0.036807235250815015, 'fraction_numerical': 0.009622462929855927, 'mean_word_length': 4.535069848661234, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Understanding Machine-learned Density Functionals

Kernel ridge regression is used to approximate the kinetic energy of non-interacting fermions in a one-dimensional box as a functional of their density. The properties of different kernels and methods of cross-validation are explored, and highly accurate energies are achieved. Accurate {\em constrained optimal densities} are found via a modified Euler-Lagrange constrained minimization of the total energy. A projected gradient descent algorithm is derived using local principal component analysis. Additionally, a sparse grid representation of the density can be used without degrading the performance of the methods. The implications for machine-learned density functional approximations are discussed.

Since the early days of quantum mechanics, it has been known that sufficiently accurate solutions of Schrödinger's equation for electrons and nuclei yield good predictions of the properties of solids and molecules [1]. But the Coulomb repulsion between electrons causes the computational cost of solving the Schrödinger equation to grow rapidly with the number of electrons, N [2]. However, as Hohenberg and Kohn proved in 1964 [3], the one-electron density may be used as the basic variable of quantum mechanics instead of the wavefunction, greatly reducing the complexity of the computational problem. This is called density functional theory (DFT) [4]. In principle, the mapping of the Schrödinger equation to one with the electron density is exact, but in practice, both the kinetic energy and the energy of the interaction between electrons must be approximated. In the original Thomas-Fermi theory [5,6], a local density functional approximation to the kinetic energy is used. However, Thomas-Fermi theory proved unsuitable for chemical and solid-state applications as it does not bind matter [7]. Shortly after the Hohenberg-Kohn theorems, Kohn and Sham (KS) [8] found a middle ground by mapping the many-body system onto a fictitious system of non-interacting electrons which reproduce the exact electron density. The main reason KS DFT became successful is because the kinetic energy of these non-interacting electrons is an excellent approximation to the many-body kinetic energy. Simple approximations to the interaction energy produce much greater accuracy and reliability compared with the standard orbital-free DFT schemes built on Thomas-Fermi theory. However, the accuracy of the results are still sensitive to the approximation of the exchange-correlation (XC) functional. In the past four decades, there has been extensive research into improving density functional XC approximations. Development of both empirical and non-empirical functionals require great intuition built on years of experience, as well as painstaking trial and error [9][10][11].
Despite the great success KS DFT has enjoyed, the computational cost scales as O(N 3 ), which is much worse than the linear scaling of orbital-free DFT [12]. Thus, there continues to be strong interest in improving upon existing orbital-free approximations to the kinetic energy [12][13][14]. A sufficiently accurate approximation to T S [n], the kinetic energy of KS electrons as a functional of the ground-state density n(r) would enable highly accurate orbital-free DFT calculations with the same accuracy as KS DFT at a fraction of the computational cost. For example, benchmark orbital-free DFT calculations are capable of treating millions of atoms in metals [15] or proteins in solvent [16]. Note that accuracy in T S beyond that of current XC approximations would be unnecessary, since all standard orbital-free DFT schemes utilize the KS decomposition of the energy, so that standard XC approximations developed for KS DFT can be utilized. However, since T S is typically comparable to the total energy of the system [4], an unavoidable problem is that a useful kinetic energy (KE) functional calls for much stricter relative accuracy than XC functionals. Additionally, accurate functional derivatives are required because one finds the ground state density by solving an Euler equation with the approximate kinetic energy functional. Continued efforts have been made in this research direction, with some notable progress [17][18][19][20][21][22][23][24][25][26][27][28]. For a review of state-of-the-art orbital-free DFT functionals, we refer the reader to Ref. [12].
In DFT, functionals typically fall into two categories. Non-empirical functionals derived from first principles tend to work well across a broad range of systems, and may exhibit systemic errors in treating certain types of interactions. Semi-empirical functionals introduce parameters that are fitted to standard data sets, and are typically more accurate with less systematic errors.
Recently, some of us applied machine learning (ML) in a completely new approach to approximating density functionals [29,30]. In a proof of principle, kernel ridge regression was used to approximate the kinetic energy of non-interacting fermions confined to a 1d box as a functional of the electron density [29]. In that work, a modified orbital-free DFT scheme was able to produce highly accurate self-consistent densities and energies that were systematically improvable with additional training data. ML algorithms are capable of learning high-dimensional patterns by non-linear interpolation between given data. These powerful methods have proved to be very successful in many applications [31], including medical diagnoses [32], stock market predictions [33], automated text categorization [34], and others. Recently, ML has been applied to quantum chemistry, including fast and accurate modeling of molecular atomization energies [35][36][37], optimizing transition state theory dividing surfaces [38], and calculating bulk crystal properties at high temperatures [39].
This new approach to density functional approxima-tion suffers none of the typical challenges found in traditional approximations, but presents many new ones. First and foremost, ML is data-driven: reference calculations are needed to build a model for the KE functional. Since every iteration in a KS DFT calculation provides an electron density and its exact non-interacting KE, reference data is relatively easy to obtain. Additionally, the ML approximation (MLA) to the KE may have thousands or millions of parameters and satisfy none of the standard exact conditions in DFT, such as positivity, scaling, and exactness for a uniform electron gas. On the other hand, the form of the MLA is completely general and thus directly approximates the functional itself, suffering none of the typical issues plagued by standard functionals starting from a local approximation. For example, some of us recently showed that an MLA for the KE has no problem accurately dissociating soft-Coulomb diatomics in 1d-a huge challenge for standard approximations [30]. However, kernel ridge regression is strictly a method of interpolation. An MLA can only be used on systems it was designed for.
In this paper, we explore the properties of the MLA derived in Ref. [29] in greater detail. In particular, we investigate the use of various kernels and their properties and the efficiency of various cross validation methods. We discuss the issue of functional derivatives of the MLA in greater detail, and explain how a modified constraint to the standard Euler equation enables highly accurate self-consistent densities, or constrained optimal densities, to be found. Additionally, a projected gradient descent algorithm is derived using local principal component analysis in order to solve the modified Euler equation. Finally, we explore the use of a sparse grid representation of the electron density and its effects on the method.

II. THEORY AND BACKGROUND
Throughout this work, we consider only noninteracting same-spin fermions in one-dimension. Thus, all electron densities n(x) are fully spin-polarized. Atomic units are used in symbolic equations, but energies are usually presented in kcal/mol.

A. Model system
Consider N non-interacting same-spin fermions subject to a smooth external potential in one-dimension, with hard walls at x = 0 and x = 1. We restrict this study to a simple class of potentials, namely a sum of 3 Gaussian dips with varying heights, widths and centers: for x ∈ [0, 1], and v(x) = ∞ elsewhere. The Hamiltonian for this system is simplyĤ =T +V , wherê T = −∂ 2 /2∂x 2 andV = v(x). We solve the Schrödinger equation for the eigenvalues j and orbitals φ j (x). As our fermions are same-spin, each orbital φ j (x) is singly-occupied. Thus, the electron density is given by and the kinetic energy is A dataset is created by randomly sampling a j ∈ [1, 10], b j ∈ [0.4, 0.6], c j ∈ [0.03, 0.1], to generate 2000 different potentials. For each potential, the system is occupied with up to 4 fermions, and the exact densities and kinetic energies are computed. Numerically, the Schrödinger equation is solved by discretizing the density on a grid: and ∆x = 1/(N G − 1) is the grid spacing. Numerov's method [40] together with a shooting method is used to solve for the eigenvalues and eigenfunctions of Eq. (2). For N G = 500, the error in our reference kinetic energies is less than 10 −7 . Fig. 1 gives a few sample densities and their corresponding potentials. The data used here is identical to that of Ref. [29]. The exact values of the parameters used in each sample is given in the supplementary information of Ref. [29]. Of the 2000 samples generated, the first half is reserved for training while the second half is reserved for testing (which we refer to as the test set).

B. Orbital-free DFT
In orbital-free DFT, T S is approximated as a functional of n(x). For our model system with non-interacting fermions, the total energy is given as for a given potential v(x). The potential is known exactly as a functional of n(x): Via the variational principle, the ground-state density is found by the Euler-Lagrange constrained search where the chemical potential µ is adjusted to produce the required particle number N . This becomes simply The density that satisfies this equation, minimizing E v [n] with the normalization constraint, is found self consistently.
Given the exact functional T [n], solving Eq. (9) will yield the exact ground-state density of the system. But in practice, T must be approximated. LetT be such an approximation, n(x) be the exact density, andñ(x) be the self-consistent density found withT . There are two measures of the error of such an approximateT [41]. The first is to compute the functional-driven error ∆T F = T [n] − T [n], which is simply the error in the KE evaluated on the exact density. The second (and much more difficult) test is to insertT into Eq. (9), solve for the approximate densityñ, and compute its error relative to the KE of the exact density . Then the density-driven error is defined as ∆E D = ∆E − ∆T F [41]. This is the additional error incurred by the approximate density. In practice, a functional which only satisfies the first test is not much use, as the ground-state density itself must also be obtained from this approximation. In orbital-free DFT, self-consistent results can be much worse than energies of KS densities, as inaccuracies in the functional derivative can cause large errors in the corresponding density. In the case of the KE functional for real systems, functional derivatives of traditional approximations can have singularities at the nuclei, making all-electron calculations very difficult, if not impossible, to converge [12]. Many of these problems can be avoided through use of pseudopotentials [12,28], but in general the solution for Eq. (9) is nontrivial.
As mentioned above, the simplest density functional approximation to T S is the local approximation [4], which for spin-polarized densities in 1d is For N = 1, the exact KE has the von Weizsäcker [17] form: As was shown in Ref. [29], the local approximation does poorly. The mean absolute error (MAE) on the test set is 217 kcal/mol, and self-consistent results are even worse at 1903 kcal/mol. A standard extension of the local approximation to a semi-local form is to add a fraction of , forming a modified gradient expansion approximation. It was shown in Ref. [29] that this did little to improve upon the local approximation.

C. Data topology and representation
Typically in ML, the data has a finite representation. For example, in Ref. [35], molecular structures are represented by a Coulomb matrix and the model predicts atomization energies. In contrast, the electronic density n(x) is a continuous function restricted to the domain [42] where H 1 (R) is a Sobolev space 1 . Although J N is infinite dimensional, in practice n(x) is expanded in a finite basis (with N G basis functions). In this work, we use a real space grid to represent the density, since our reference calculations are done using the same grid. We use the L 2 inner product and norm between densities n i (x), n j (x) (In actual calculations, all densities are represented on a finite basis, and thus will have have a finite L 2 -norm).
Since the ML algorithm is expressed in terms of this inner product, the results are independent of the specific representation used as long as the basis is converged. Even with a truncated basis, J N is still highdimensional and applying ML to learn the KE of all densities in J N would not be feasible. Fortunately, we are only interested in a subspace of J N related to a specific class of potentials (e.g. Gaussian dips), which greatly reduces the variety of possible densities. In general, let the potential v(x) be parametrized by the parameters We define the density manifold M N ⊂ J N as the set of all densities that come from these potentials with a given particle number N . In general, M N is a d-dimensional manifold. The training densities, n j (x) for j = 1, . . . , N T , are sampled from M N . In the present work, the external potential has 9 parameters, and thus d is at most 9.

D. The kernel trick and feature space
In finding the structure of low-dimensional data, it is often sufficient to optimize parametrized non-linear forms (e.g., using a polynomial to fit a sinusoid). For high-dimensional, nonlinear data this becomes increasingly difficult. In kernel-based machine learning, the approach is to transform the data itself non-linearly to a high-dimensional space known as feature space, such that the data becomes linear [31,[43][44][45][46]. Fig. 2 illustrates data points that lie on a circle in the Cartesian plane. As shown, the data becomes linear on transformation to polar coordinates, and linear regression can subsequently be used to fit the data. Transforming back to Cartesian coordinates recovers the non-linearity. Let the data points belong to a vector space χ, also called input space, and let Φ : χ → F be the map to feature space F . Assuming we wish to apply a linear method such as regression in feature space F , we note that regression can be expressed solely in terms of the inner product between feature vectors Φ(x) and Φ(y), where x, y ∈ χ. We define the kernel k such that The kernel can generally be thought of a measure of similarity between data, but must satisfy Mercer's condition: Mercer's theorem [47] guarantees the existence of a feature space F , which is a reproducing kernel Hilbert space [48]. Since the linear algorithm in F may be expressed in terms of the kernel in Eq. (14), Φ need never be explicitly computed. This procedure is known as the kernel trick, and enables easy nonlinearization of all linear scalar productbased methods that can be expressed via an inner product [49].

E. Kernel ridge regression
Kernel ridge regression is a nonlinear version of regression with a regularization term to prevent overfitting [50]. Our MLA for the KE has the form where N T is the number of training densities, α j are weights to be determined, n j are training densities and k[n, n j ] is the kernel. The weights are found by minimizing the quadratic cost plus regularization , and λ is called the regularization strength. The second term penalizes weights with large magnitudes in order to prevent overfitting. 2 By setting the gradient of Eq.17 to zero, minimizing C(α) gives where I is the identity matrix and . The hyperparameters, which include the regularization strength λ and the parameters of the kernel such as the length scale σ, are found via cross validation (see [36] and Sect. III).
The choice of the kernel will depend on the given data. Some kernels are designed to be generally robust and applicable (e.g., the Gaussian kernel), while others are designed for a specific type of data (see e.g. [31,52,53]). A good choice of kernel can reflect the characteristics of the data (see [54]). In Ref. [29], we chose the Gaussian kernel 2 The regularization term accounts for the possibility of noisy data (e.g. experimental data), and imposes certain smoothness conditions on the model (see [51]). Our reference data is deterministic and thus noise-free in this sense, but, because the precision of our calculations is limited, we may consider the numerical uncertainty to be noise. where σ is the length scale. Since the density is represented on a uniform grid, the L 2 -norm can be approximated by 3 where x l is given by the grid defined in Eq. (5). This approximation becomes exact as ∆x → 0. Fig. 3 shows the range and distribution of Euclidean distances between all pairs of densities and KE of all densities in the dataset with N = 1. Ordinary linear regression models frequently employ a bias term to account for the fact that the data might lie away from the origin. Without this term, the regression line is forced to go through the origin, causing a systematic error if the data does not. The bias term can be implemented directly, or by centering the samples such that the mean is zero. Fig. 4 illustrates the transformation to feature space for an example 1d data set and linear regression in feature space. If the data is centered in feature space, the bias term is unnecessary.
Here, we center the densities in features space such that For simplicity, all equations given in this work assume that the data is centered (i.e. k =k). In fact, kernels such as the Gaussian kernel Eq.(19) whose induced reproducing kernel Hilbert space on a bounded domain is dense in the space of continuous functions on this domain do not require centering [55].

A. Kernels
Model selection refers to the process of selecting a kernel and the corresponding hyperparameters. In kernel ridge regression, this includes the regularization strength λ and the kernel parameters (e.g. in the Gaussian kernel, the length scale σ). Table I lists some standard kernels.
Radial basis function (RBF) kernels, which include the Gaussian, Cauchy, and Laplacian kernels, all behave similarly and tend to work for a broad range of  problems. Other kernels work well for specific data structures [31,52,53] and regularization properties [49].
Fig . 5 shows the contours of the functional-driven MAE over the test set as a function of the regularization strength λ and the kernel parameter σ. We see that the qualitative behavior is similar for the Gaussian, Cauchy and Laplacian kernels. In the left region (where the contour lines are vertical), the length scale σ is much smaller than the distance between neighboring training densities. Thus the RBF-type kernel functions centered at each training density have minimal overlap, yielding a poor approximation to the KE functional. The kernel matrix becomes nearly unity, and the regularization λ has negligible effect. On the right side of the contour plot, the length scale is comparable to the global scale of the data. In these regions, the kernel functions are slowly varying and do not have enough flexibility to fit the nonlinearity in the data. The region with minimum MAE lies in the middle. The Gaussian and Cauchy kernels both give the same performance, with errors less than 1 kcal/mol in the middle region (enclosed by the dashed line), while the Laplacian kernel behaves poorly in comparison. This is likely due to the cusp in the form of the kernel, which cannot fit the smooth KE functional.

B. Optimization of hyperparameters
After picking a kernel family, the values of the hyperparameters must be chosen. Ideally, we select the hyperparameters such that the generalization error, which is the error not only on our training set but also on all future data, is minimal. The out-of-sample error must be estimated without looking at the test set (the test set is never touched during model selection, so that it can give a true test of the final performance of the model) [31,36]. This procedure, known as cross-validation, is essential for model selection in preventing overoptimistic performance estimates (overfitting).
Various schemes for cross validation exist [31,36,56], but all obey a basic principle: the available data is subdivided into three parts: the training, validation and the test sets. The ML model is built from the training set and the hyperparameters are optimized by minimizing the error on the validation set (Fig. 6). The test set is never touched until the weights and hyperparameters have been determined. Then and only then, the generalization ability of the model can be assessed with the test data ( [31,36], see also [57]). Typically, the data is shuffled to ensure its random distribution between training and validation division. This can be repeated with different subdivisions. A few schemes, which will be analyzed for our kinetic energy functional estimation problem, are described below. For each scheme, a test set of 1000 samples is used to estimate the generalization error after the ML model is selected.
Simple cross validation: The training data (N T samples) is randomly divided into a training set of 70% and a validation set of 30%. The hyperparameters are optimized by minimizing the MAE on the validation set.
k-fold cross validation: Step 1: The N T training data is randomly divided into k bins. Step 2: The jth bin is used as the validation set and the remaining k − 1 bins as training set. The model is built on the training set and the hyperparameters are selected by minimizing the MAE on the validation set.
Step 3: Repeat step 2 k times such that all bins have been used as validation sets. We will then have k models in total and the final hyperparameters are selected as the median over all models.
Because the mean cross validation error still depends on the initial random partitioning of data in cross validation, we repeat the procedure with different subdivisions [31].
Leave-one-out: Leave-one-out (LOO) is a special case of k-fold cross validation, when k = N T . Thus each bin contains only one sample.
Typically, it is better to leave out as little data as possible to exploit the statistical power in the data. Simple cross validation is computationally expedient, but wasteful since not all training data participates in the optimization. k-fold cross validations are used in situations where data is very limited, or expensive to collect. Leaveone-out is often used with limited data and it becomes computationally intensive if N T is large. k-fold cross validation gives a good balance on all counts.
For the Gaussian kernel, Fig. 7 shows the MAE on the test set with the hyperparameters optimized with different cross validation methods. With 120 training densities, all schemes give a similar MAE, despite the large variations in σ and λ. This means that multiple models exist that give comparable performance. As expected, the only variations in MAE occur for more limited data. Fig. 5 shows how 10-fold cross validation performs in selecting hyperparameters that generalize well to the test set, for a few kernels. The gray dots represent the optimal parameter choice for each repetition, and the black dot is the median over all repetitions. In this case, the global  Table II shows the optimal hyperparameters and functional driven errors for the kernels listed in Table I. Some optimum values for the Gaussian kernel are listed in Table III. Detailed information with optimum values of other kernels are shown in supplementary material.

IV. RESULTS AND DISCUSSION
In the main work of this paper, we test in greater detail some of the methods that were introduced in Ref. [29] using only the Gaussian kernel, as it performs the best.

A. Errors on exact densities
In Table III, we evaluate our MLA, constructed using the first N T training densities in our data set, on the exact densities of the test set and compute the errors . The Gaussian and Cauchy kernels give the best performance. For the Gaussian kernel with N = 1 chemical accuracy is achieved (i.e. MAE less than 1 kcal/mol) at N T = 60. Just as we saw in Ref. [29], the performance is systematically improvable with increasing number of training densities. The Laplacian kernel gives a mean absolute error (MAE) of 6.9 kcal/mol at N T = 100 (still better than LDA), which improves as N T increases. On the other hand, the performance of the wave kernel does not improve as N T increases (see supplemental information). This indicates the form of the wave kernel is not flexible enough to fit the form of the KE functional.

B. Sparse grid
Note that the choice of N G used in the reference calculations is needed to converge our reference energies and densities, but may be larger then the grid needed to "converge" our ML functional. As the ML model depends only on the inner product between densities, this will typically converge much faster than, e.g. Numerov's method. To demonstrate this, we define a "sparse" grid, {x s(j−1)+1 |j = 1, . . . , N G /s}, using every sth point in the grid (we only choose s such that N G is divisible by s). Fig. 8 shows that performance of the model is unaffected until N G is reduced to about 10 grid points. The model is cross-validated each time, but the hyperparameters change only slightly. Thus, ML can accurately learn the KE functional with a far less complete basis than is required to accurately solve the Schrödinger equation. This is possible because we have restricted the learning problem to a simple type of potential with a limited range of possible densities and energies. The underlying dimensionality of the data is about 9, comparable to the number of parameters that determine the potential. The model needs only enough degrees of freedom in the representation of the density to distinguish between densities, but no more. Thus it is no coincidence that the minimum grid required is comparable to the dimensionality of the data (i.e. the dimensionality of the density manifold M N ).
However, we also need a sufficiently fine grid to compute the integral in Eq. (7) to the desired accuracy. In the problem shown here, the dimensionality of the data is relatively small, and will increase for larger systems (e.g. real molecules with many degrees of freedom). In general, however, we need to consider both factors in choosing a suitable basis. But, we may be able to use a basis that is more sparse than that of the reference data, which would greatly reduce the computational cost of the method.

C. Challenge of finding density
Thus far, we have focused on the discussion of the performance of the MLA evaluated on exact densities (i.e. the functional-driven errors). However, in order for a functional to be useful, it must also predict the groundstate density. As discussed previously, an accurate functional derivative is necessary in order to solve Eq. (9) and yield an accurate density. The functional derivative of our MLA is given by: where, for the Gaussian kernel, δk[n, n j ]/δn(x) = (n j (x) − n(x))k[n, n j ]/σ 2 .
In Fig. 10, we plot the functional derivative of our model compared with the exact derivative. The model displays   gives the error of KE evaluated on constrained optimal densities in kcal/mol and (c) the corresponding re-cross validated hyperparameters λ and σ. The MAE is completely unaffected as NG is reduced until approximately NG = 10, when it jumps sharply.
a highly inaccurate functional derivative, with a huge amount of apparent "noise", as was found in Ref. [29].
What is the source of this noise? In general, if the underlying dimensionality of the data is much less than the dimensionality of J N (which in this case is essentially infinite), ML will be unable to capture the functional derivative. The functional derivative contains information on how the KE changes along any direction, but ML cannot learn this because it only has information in directions in which it has data (i.e. along M N ). Fig. 11 illustrates the problem: standard minimization techniques will rapidly exit the "interpolation" region in which the MLA is expected to be accurate. The MLA is only given information about how the KE changes along the density manifold M N . In the many dimensions orthogonal to M N , the MLA produces an inaccurate derivative (each of these dimensions produces a large relative error since no data exists in these directions; the sum over many dimensions creates a large total error in functional derivative). A standard gradient descent will quickly venture off of M N into regions of J N where the model is guaranteed to fail. Fig. 9 shows the deviation of self-consistent density if the search is not constrained to M N . To fix this, we further constrain the minimization in Eq. (9) to stay on M N . The Euler-Lagrange minimization for the ground-state density can be expressed as where g is any function that is zero on M N and positive elsewhere. the minimizing density of this equation the constrained optimal density. It may be solved self-consistently in the same sense of solving the standard Euler equation. However, the g[n] which exactly gives the density manifold is unknown. In the next section, we develop an approximation which attempts to reconstruct the density manifold from the training densities.

D. Manifold reconstruction using principal component analysis
Our aim is to reconstruct M N locally around a given density n(x), which is assumed to be on the density manifold. A simple approach is to approximate M N as locally linear, using principal component analysis (PCA) to determine the tangent space empirically from the training densities. This will work as long as there are enough training densities covering the density manifold. First, we define a weighted average density around density n: This generalized average is weighted by the function ω( n − n ) that only depends on the distance from n(x) to n (x), ω j = ω( n − n j ), and Ω = N T j=1 ω j . Note that n (x) refers to the density n evaluated at x and not the derivative of n with respect to x.
The locality of the method comes from the choice of ω. For standard PCA, the choice is ω(r) = θ(R − r), where θ is the Heaviside function, and R is the distance from n to the m-th nearest training density. This equally weights the nearest m training densities, and ignores all other training densities. This choice was used in Ref. [29]. Here, we choose a slightly smoother weighting function: Next, PCA is performed by spectral analysis of the empirical covariance operator [58], based on the weighted average value around n(x). We define the centered neighborhood byñ j (x) = n j (x) −n(x). In this problem, densities are represented on a grid with N G = 500 points, so let n = (n(x 1 ), . . . , n(x N G )) be the vector representation of n(x). The covariance matrix Γ ∈ R N G ×N G is with eigendecomposition The eigenvalues are ordered such that λ j > λ j+1 . The eigenvectors u j are called principal components (PCs), and give the directions of maximum variance in the data. We define the variance lost in keeping d PCs as η = 1 − d j=1 λ j N G j=1 λ j . In this case, there is little to no variance in directions orthogonal to the tangent space of M N , and maximum variance in directions aligned with the tangent space. Thus, the first d PCs form a basis for the tangent space, where d is the dimensionality of the density manifold (and tangent space). The projection operator onto this basis is: The tangent space using PCA is given by Finally, we choose the PCA approximation to the constraint g[n] in Eq. (24) as the squared distance from n to tangent plane T PCA [n]: The PCA approximate density manifold M PCA is then defined implicitly by g PCA [n] = 0. The process is illustrated in Fig. 12. In the next section we develop a projected gradient descent method to solve Eq. (24).

E. Projected gradient descent algorithm
For a given ML approximation to the KE functional, the algorithm to minimize the functional in Eq. (24) to find a constrained optimal density is as follows (see Fig.  13). Choose an initial guess for the density, n 0 ∈ M N (e.g., a training density): 1. Evaluate the functional derivative at n = n t .
3. Project the functional derivative onto the tangent space (see Fig. 13), and take a step: where is a constant such that 0 < ≤ 1. If convergence is unstable, reduce , trading stability for speed of convergence.
4. To ensure the constraint remains satisfied, we subtract the (weighted) mean of the training densities in the local neighborhood: We iterate these steps until convergence is achieved. We measure convergence by setting a maximum iteration step and tolerance threshold. If the total energy difference is smaller than tolerance within max iteration step, the density is converged. If no solution is found, is reduced.  13. Schematic of the projected gradient descent. The functional derivative is projected onto the tangent space of the data manifold MN at nt (dashed line). Next, a step is taken along the projected functional derivative to n t in the direction of lower energy. Finally, g[n] is minimized orthogonal to the tangent space to ensure the minimization stays on MN .

F. Errors on constrained optimal densities
With this new constrained minimization procedure via a projected gradient descent, we solve for the constrained optimal density for each test sample. We report the errors in the total energy and KE relative to the exact density in Table III. In general, we expect these errors to be worse on the MLA evaluated on exact densities-by roughly a factor of 10. However, errors on constrained optimal densities decrease at the same rate with more training data, so an accuracy of 1 kcal/mol in KE is achieved with 150 training samples for N = 1, now on constrained optimal densities. Additionally, errors are of similar magnitude for multiple particles. In the last row of Table III, we combine the training data from each N (100 training densities per N value) into one model. This combined MLA gives roughly the same error as each individual model. This is because, due to the locality of the Gaussian kernel, the training densities from each N are well separated (orthogonal in feature space) and the individual models are unaffected.
In the projected gradient descent, there are two PCA parameters that must be chosen: m, the number of nearest neighbors and d, the number of PCs to form the projection operator. Fig. 14 shows the MAE evaluated by the constrained optimal density and variance lost as a function of the number of PCs d with m = 20. The MAE decreases initially as d increases as more PCs capture the local structure of the density manifold. As can be seen, d = 4 or 5 gives an optimal reconstruction of the tangent space of the manifold. As d increases further, the noise that was removed is re-introduced into the projection, causing the gradient descent algorithm to fail. For d = 7, many of the constrained searches do not converge. Table IV reports the errors of the model evaluated on constrained optimal densities for N T = 40 and N T = 100, giving a rough optimization of the PCA parameters. Although the potential which generates M N has 9 parameters in this case, we observe that the optimal choice of d is only 4. This is because the data used The PCA projection will compensate by removing those directions. Thus, the effectiveness of our method depends on the sampling on the manifold.

V. CONCLUSION
In this work, we have explored in much greater detail the methods presented in Ref. [29], in which ML methods were used to directly approximate the KE of a quantum system as a functional of the electron density, and used this functional in a modified orbital-free DFT to obtain highly accurate self-consistent densities and energies.
We used a simple model as a proof of principle, to investigate how standard methods from ML can be applied to DFT. In particular, we tested a variety of standard kernels used in ML, and have found that the Gaussian kernel gives the lowest errors (the Cauchy kernel also achieves similar performance). All cross validation schemes that were tested gave similar predictions of hyperparameters that achieved low generalization error on the test set. Our results highlight the importance of an appropriate choice of kernel, as some of the kernels tested gave strikingly bad performance. With the construction of the L 2 norm that was used in the kernels, the method is basis set independent (as long as a complete basis is used). However, the ML method is capable of learning accurate KEs using a sparse grid (i.e., an incomplete basis). Using a sparse representation for the density without losing accuracy would speed up calculations further. These re-sults warrant further exploration and will be the subject of future work.
We explained the origin of the noise in the functional derivative and developed a constrained search over the density manifold via a modified Euler equation, effectively projecting out the noise. We also introduced a local approximation to the manifold using PCA, and solved for constrained optimal densities using a projected gradient descent algorithm. This worked well for our prototype system, yielding highly accurate constrained optimal energies and densities.