The Levy-Lieb embedding of density functional theory and its Quantum Kernel: Illustration for the Hubbard Dimer using near-term quantum algorithms

Density Functional Theory (DFT)  [1, 2] is presently the dominant paradigm for material-specific electronic structure simulations in computational materials science [3]. Originally proposed by Hohenberg and Kohn (HK) [1], DFT establishes the one-body density as the basic variable in interacting many-body problems with fixed interaction strength subject to time-independent spatially local external potentials [4]. The electronic Hamiltonian in the Born-Oppenheimer approximation is the most prominent example of such a setup. Given a problem in this class, all ground state observables can in principle be represented as explicit functionals of the ground state one-body density. However, while the density is informationally adequate, deriving explicit density functionals that are universally accurate is a hard problem in the same complexity class as computing generic ground states [5]. The practical success of DFT instead stems from the development of extremely low cost yet usefully accurate [6] approximate models for the quantum mechanical correlation energy in the form of highly compact density functionals [7] circumventing the need to represent many-electron wavefunctions explicitly on classical computers at simulation time. Benefits of the compactness of a density functional representation can also extend to quantum systems where efficient classical heuristics for the wavefunction exist [8, 9]. In the modern context the encodings implied by DFT may be viewed as simply the most condensed description within a hierarchy of functional theories based on reduced density matrices [5, 10].

With the advent of quantum computers the prospect of representing many-electron states on quantum hardware has important implications for DFT and other functional theories [5, 11, 12, 13, 14]. The classical intractability of interacting many-electron wavefunctions partly provided the motivation for functional theories [15] but at the same time also constrained the development of approximate functionals [9]. For instance the study of exact forms of DFT has to date been largely restricted to small model systems where systematic comparisons with wavefunction calculations are possible [9] and pure density functionals for a wider variety of observables are known in model systems than are available for deployment in realistic simulations of matter [9, 8, 16, 17]. Furthermore, the relationship between approximate functionals and underlying many-electron states is often obscured in practice even though the mapping can be recovered ex post facto at significant computational cost [18]. Finally and most relevant to this work is the fact that outside of a few pioneering efforts [19, 20, 21, 22, 23, 24] the explicit density to wavefunction mappings implied by DFT [25, 26, 27, 28, 29] are not widely pursued in numerical studies on classical computers and one typically encounters density functionals only after they are approximated either using formal, empirical or machine learning methods [7, 30]. The availability of quantum devices may alleviate this constraint in future. Therefore in this work we calculate the density to wavefunction mapping [28, 29, 31, 19, 20, 21, 24], as formulated by Levy [25, 26] and Lieb [27] using a density-constrained variational quantum eigensolver (VQE) [32, 33] scheme and demonstrate density variational minimization to identify the ground state of a paradigmatic fermion lattice problem, the Hubbard dimer [16, 17]. Furthermore, we make use of the density-wavefunction map to define a quantum kernel [34, 35, 36, 37] which we call the Levy-Lieb quantum kernel and use it to learn observable functionals in such a way that the relationship to the underlying state space is readily apparent. Our work contrasts with previous research efforts which explored proposals for the popular Kohn-Sham DFT [2] in both near-term and fault-tolerant settings [11, 12, 13, 14]. The rest of the article is organized as follows. In section II, we introduce the Levy-Lieb procedure which defines the density to wavefunction map and in section III we outline a concrete VQE based implementation of the same for lattice systems. Section  IV discusses results of an explicit density variational search for the ground state of the Hubbard dimer. In section V we introduce a fidelity based Levy-Lieb quantum kernel (LLQK) and illustrate its use in machine learning density functionals. We summarize our conclusions in section VI.

Consider an interacting $N$-particle system described by a Hamiltonian operator of the form

where $\hat{T},\hat{W}$ are the kinetic and inter-particle interaction operators and the external potential $\hat{v}$ is a local one-body operator. Following Eschrig [31], we define the set of $N$-representable densities as

with $\bm{L}^{2}$ denoting the space of square-integrable functions, and the set of $N$-particle wavefunctions as

The Levy-Lieb density functional $F_{\mathrm{LL}}[n]$ is then defined by

i.e., given a density $n\in\mathcal{J}_{N}$, the Levy-Lieb procedure involves searching for the infimum of $\langle\Psi|\hat{T}+\hat{W}|\Psi\rangle$ over the restricted subset of wavefunctions $\Psi\in\mathcal{W}_{N}$ constrained to yield the chosen density $n$. We use the notation $\Psi\rightsquigarrow n$ to indicate $\Psi$ yields $n$ and define $\mathcal{W}^{n}_{N}\equiv\{\Psi\in\mathcal{W}_{N}|\Psi\rightsquigarrow n\}$. Note that the above procedure does not involve the external potential $\hat{v}$ and at the end of the search process we have access to both the value of $F_{\mathrm{LL}}[n]$ and a $\Psi\rightsquigarrow n$. For general $n\in\mathcal{J}_{N}$ the obtained $\Psi$ belongs to a sub-manifold of states on $\mathcal{W}^{n}_{N}$ with the same expectation $\langle\Psi|\hat{T}+\hat{W}|\Psi\rangle$ [31] but in the absence of such degeneracies on $\mathcal{W}^{n}_{N}$, $n$ uniquely determines $\Psi$ [39]. In particular, on the subset of densities that correspond to non-degenerate ground states, the mapping is bijective due to the Hohenberg-Kohn theorem [1, 31, 39]. In all cases, $F_{\mathrm{LL}}[n]$ is well-defined. The density variational principle follows directly from the Levy-Lieb procedure. Consider the set $\mathcal{W}_{N}$ and define the equivalence relation

This leads to a partitioning of $\mathcal{W}_{N}$ into disjoint subsets $\mathcal{W}^{n}_{N}$ labelled by the density $n$ as $\Psi_{i}\in\mathcal{W}^{n}_{N}~{}\mathrm{iff}~{}\Psi_{i}\rightsquigarrow n$. On a particular $\mathcal{W}^{n}_{N}$, consider minimizing the expectation of the Hamiltonian from equation 1 with external potential $\hat{v}$

Since $\mathcal{W}_{N}=\bigcup\limits_{n\in\mathcal{J}_{N}}\mathcal{W}^{n}_{N}$, for the ground state energy we have

Thus the ground state energy $E[\hat{v}]$ can be obtained by minimizing $E[\hat{v},n]=F_{\mathrm{LL}}[n]+\int dx~{}v(x)n(x)$ over $n\in\mathcal{J}_{N}$. It follows further that if $n_{\mathrm{GS}}$ is the ground-state density then through equation 4 on $\mathcal{W}^{n_{\mathrm{GS}}}_{N}$, it determines the ground-state manifold $\Psi_{\mathrm{GS}}$. We note that the Hohenberg-Kohn (HK) functional is a restriction of $F_{\mathrm{LL}}[n]$ to the space of ground state densities [31].

As formulated by Cioslowski [19, 20, 21] and previously illustrated by Mori-Sanchez et al [24] on classical computers, the minimization procedure of equation 7 defines an exact numerical density functional calculation for the ground state provided $F_{\mathrm{LL}}[n]$ and it’s derivative ${\delta}F_{\mathrm{LL}}[n]/{\delta}n$ are calculated explicitly through constrained search. In this work we calculate these quantities for discrete lattices by employing parameterized quantum circuits executed on a statevector simulator. In the following sections when discussing formal statements unrelated to specific implementation details we will indicate functionals as $A[n]$ using square brackets and italicized $n$ as above but when referring to discretized implementations where parametric differentiation is substituted for formal functional differentiation [19, 20, 21, 40, 24], we will indicate functional dependence via $A(\mathbf{n})$ using parenthesis and bold $\mathbf{n}$ notation for the density.

Since first being proposed in 2014, the variational quantum eigensolver (VQE) [32] has been widely adopted on near-term quantum computers as a versatile hybrid classical-quantum algorithm for optimization tasks employing parameterized ansatze. The incorporation of constraints into VQE variational optimization has also been discussed previously [41, 42]. For a recent comprehensive review of VQE we refer readers to reference [33]. For our purpose, since equation 4 implies a search over wavefunction space restricted to a specified density sector, we augment the usual VQE with a density-constraint. If the Hamiltonian takes a form where the potential term appears as a local one body operator $\hat{v}=\sum_{i}v_{i}\sum_{\sigma}\hat{c}^{\dagger}_{i\sigma}\hat{c}_{i\sigma}$ in second-quantization, then the one-body density operator of interest needed to specify the constraint is $\hat{n}_{i\sigma}=\hat{c}^{\dagger}_{i\sigma}\hat{c}_{i\sigma}$ which is easily computed on any discrete index set $\{i\sigma\}$ enumerating basis modes $i$ with spin index $\sigma$. In the literature, this choice of relevant potential and density on a lattice goes by the name site-occupation functional theory [9] but the same form also appears in connection with discretized real space grids [24]. Then for $N$ electrons on $M$ sites and parameterized ansatz states $|\psi(\theta)\rangle$ with real parameters $\theta$, the Levy-Lieb procedure consists of minimizing the expectation

subject to the constraints

where $\{n_{i},i=1..M\}$ are prescribed site occupations and additionally we have

The steps involved in density-constrained VQE (DC-VQE) are outlined in table 1. We expect based on the work of D’Amico et al [23] that for $v$-representable densites, i.e., densities that correspond to $N$-particle ground states in some potential $v$, small errors in computing densities do not lead to large errors in the wavefunction as nearby densities get mapped to nearby wavefunctions in metric space and densities on lattices are $v$-representable under very reasonable assumptions [43, 44]. Thus the constrained minimization process should be robust to small numerical errors in computing densities.

Once a subroutine for calculating $F_{\mathrm{LL}}(\mathbf{n})$ is available, density variational minimization for a given external potential can be implemented as outlined in table 2 wherein, following earlier implementations for classical computers [19, 20, 21, 24], we compute the derivative ${\delta}F_{\mathrm{LL}}(\mathbf{n})/{\delta}\mathbf{n}$ using finite differences while maintaining normalization of the density. The next section outlines a specific implementation of DC-VQE and density variational search for the case of the asymmetric Hubbard dimer.

The asymmetric Hubbard dimer has been explored extensively as a paradigmatic model system within lattice density functional theory [16] and time-dependent density functional theory [17] and serves as an ideal test case to explore near-term variational quantum algorithms and quantum machine learning in the context of DFT. The relevant fermion Hamiltonian

features a local multiplicative external potential term as the source of inhomogeneity in the model. It is further characterized by two dimensionless parameters $U/t$ and $\Delta v/t$ where $\Delta v=v_{1}-v_{2}$ is the potential asymmetry between the two dimer sites. To enable direct comparison with previous benchmark results [16, 17] we fix $2t=1$ and vary $U$, $\Delta v$ across different calculations. Additionally, since for the on-site occupations $\mathbf{n}=[n_{1},n_{2}]$ we have $n_{1}+n_{2}=N$, a single parameter $n_{1}$ is chosen to specify the density distribution. We consider the case of $N=2$ particles in four site-localized spin orbitals and restrict our attention to the lowest energy spin-singlet ($S^{2}=0,S_{z}=0$) which is sufficient for analyzing ground state DFT of the dimer [16]. The above setup leads via the Jordan-Wigner mapping [45] to a 4-qubit problem on a quantum computer. For this simple model we find that a generalized unitary coupled cluster [46] (UCC) wavefunction ansatz featuring single and double excitations along with two Trotter steps is sufficient to reproduce the exact-diagonalization result for the ground state energies over the parameter range of interest in conjunction with VQE [32, 33] (See fig 1). For a general discussion of specialized ansatze for lattice models in different dimensions we refer the readers to recent literature [33]. In our study the UCC circuit ansatz which features ten parameters as shown in figure 1(a) is set up using the Qiskit [38] framework and executed on a statevector simulator backend. VQE classical parameter optimization is performed using the L-BFGS-B [47] optimizer as implemented in the SciPy [48] package. As shown in fig. 1(b) for a range of $\left(U,\Delta v\right)$ parameters, energies of the dimer spin-singlet ground-state as obtained from VQE are within $10^{-6}$ of corresponding exact-diagonalization values and in good correspondence with figure $3$ of reference [16]. We then proceed to employ the same circuit ansatz to prepare states for evaluation within density-constrained VQE simulations of the dimer.

As outlined in section III and table 1, density-constrained VQE (DC-VQE) is intended as a subroutine to calculate the Levy-Lieb functional $F_{\mathrm{LL}}(\mathbf{n})$ for a specified density $\mathbf{n}$ and makes no reference to the external potentials $\mathbf{v}$. Thus with the operator sum $\hat{T}+\hat{W}$ and a set of site-occupation difference operators $\mathbf{D}$ as defined in table 1 specified, the UCC ansatz parameters are optimized with VQE under an occupation-preserving constraint (equation 9) to yield $F_{\mathrm{LL}}(\mathbf{n})$. We implement the constraint in practice by setting up the SLSQP [49] constrained optimizer from SciPy [48] to request site-occupation measurements from the quantum simulator. For the dimer case, the DC-VQE simulated result for $F_{\mathrm{LL}}$ which is just a function of $n_{1}$ is plotted in figure 2 for different values of $U$. As noted previously by Carrascal et al [16], there is no known analytical expression for $F_{\mathrm{LL}}(n_{1})$ but it is a monotonic function on the interval $n_{1}\in(0,1]$ and approaches a bound $U$ as $n_{1}\rightarrow 0$. Furthermore since $\frac{dF_{\mathrm{LL}}(n_{1})}{dn_{1}}=\frac{\Delta v}{2}$ [16] its slope increases sharply as $n_{1}\rightarrow 0$ but remains finite for physical potentials. It is apparent that the DC-VQE result for $F_{\mathrm{LL}}(n_{1})$ exhibits the expected features.

In connection with DC-VQE it is instructive to analyze the degree to which the occupation-preserving constraint of equation 9 is respected during the optimization process for $F_{\mathrm{LL}}(n_{1})$. This is shown in figure 3(a) where we plot the constraint vector magnitude $||\mathbf{d}||$ as defined in table 1 along a typical optimization trajectory. We see that starting from a small value of $||\mathbf{d}||$ as the initial state is prepared to lie near $\mathcal{W}^{\mathbf{n}}_{N}$ the SLSQP optimizer maintains a low baseline value of $\sim 10^{-7}$ for $||\mathbf{d}||$ but with isolated spikes in between where it assumes larger values signifying departures of the trial wavefunction from $\mathcal{W}^{\mathbf{n}}_{N}$ at some instances. Therefore our approach for finding $F_{\mathrm{LL}}(n_{1})$ is not a strict constrained-search in the sense of Levy [25, 24] but since $||\mathbf{d}||$ is always small at convergence we ensure that the optimal state lies numerically close to $\mathcal{W}^{\mathbf{n}}_{N}$ as required.

With a means of calculating $F_{\mathrm{LL}}(n_{1})$ at hand, we then proceed to perform density functional simulations for a given external potential asymmetry $\Delta v$ to identify the dimer ground state. This involves a variational search over occupation-numbers as outlined in table 2 combined with on the fly calculation of $F_{\mathrm{LL}}(n_{1})$ using DC-VQE and it’s derivative $F_{\mathrm{LL}}^{\prime}(n_{1})$ using finite differences. Thus the density variational minimizer repeatedly calls DC-VQE at occupations $\mathbf{n}$ it encounters along the optimization trajectory. To initialize the DFT simulation we solve the non-interacting problem with $U=0$ and use the resulting occupation numbers as the starting guess. One could also for instance employ an approximate exchange-correlation potential $v_{xc}$ [9, 16] in the dimer Hamiltonian to provide an initial density guess that will in most cases be closer to the true ground state density than the non-interacting guess. For the given $\Delta v$ we also compute the regular VQE ground state energy and density as the benchmark result. In figure 3(b,c) we show convergence with respect to VQE of the DFT energy and site occupation $n_{1}$ as a function of the number of $F_{\mathrm{LL}}(n_{1})$ evaluations which represent the most intensive part of the simulation. We set $\Delta v=1$ and show trajectories for $U=1$ and $U=5$. In this parameter regime the correlation energy is typically sizable (see figure 6 of reference [16]) while at the same time the non-interacting and exact densities are expected to differ significantly. We find that around 20-30 evaluations of $F_{\mathrm{LL}}(n_{1})$ are needed to converge the DFT energy and occupations to within $10^{-4}$ and $10^{-3}$ of the respective VQE reference values.

To conclude this section we note that while Levy-Lieb constrained-search [25, 26, 27, 28, 29] clarifies the coarse-graining step involved in switching our description of quantum systems from wavefunctions to densities, DFT based on repeated real-time evaluations of the exact $F_{\mathrm{LL}}[\mathbf{n}]$ using explicit constraints would be inefficient relative to unconstrained minimization in wavefunction space. Furthermore constrained optimization of parameterized quantum circuits is generally non-convex [50, 33, 36, 37] and gradient based optimizers may get stuck in one of many local minima and therefore one needs to employ additional heuristics to locate the global minimum of the quantity being optimized. In the dimer case investigated here, similar to the case of variational quantum deflation on small molecules as noted previously [51], randomly initializing different starting guesses for the UCC ansatz parameters proved sufficient to identify the correct optimum for $F_{\mathrm{LL}}(\mathbf{n})$. For larger problem sizes more sophisticated global optimization schemes will be needed. Ideally specialized density-preserving parameterized ansatze that efficiently explore a specified density sector of Hilbert space would be desirable.

Restricting ourselves to non-degenerate settings where $n$ uniquely determines $\Psi$ [39] as outlined in section II, we can think of the Levy-Lieb mapping as implementing a particular feature map [35, 34, 37, 36] encoding densities into the feature space $\mathcal{F}_{\mathrm{LL}}$ of pure-state density matrices $\rho_{\mathrm{LL}}[n]\equiv|\Psi_{\mathrm{LL}}[n]\rangle\langle\Psi_{\mathrm{LL}}[n]|\in\mathcal{F}_{\mathrm{LL}}$. Additionally for the same density $n$, different embeddings can be realized through different choices for the operator sum $\hat{T}+\hat{W}$. So for a specified $\hat{T}+\hat{W}$, this allows us to define a fidelity based Levy-Lieb quantum kernel (LLQK)

for densities $n,n^{\prime}\in\mathcal{X}_{N}\subset\mathcal{J}_{N}$ where $\mathcal{X}_{N}\equiv\{n\in\mathcal{J}_{N}|n\iff\Psi_{\mathrm{LL}}\in\mathcal{W}^{n}_{N}\}$ is the set on which the Levy-Lieb map is one-to-one. Note that this restriction is not needed if one only wishes to learn the functional $F_{\mathrm{LL}}[\mathbf{n}]$ alone and furthermore in situations with degeneracies one may optionally provide additional quantum numbers to pick a particular state on the degenerate manifold to uniquely specify the LLQK. We can use the LLQK in a machine learning model of the form

where the goal is to learn the optimal operator measurement $\mathcal{O}$ or equivalently the coefficients $\alpha_{m}\in\mathbb{R}$ given a set of $M$ labelled training data $\{n_{m},f_{\mathcal{O}}[n_{m}]\}$. In particular since ground state density functionals are generated as expectation values of the form

where $\hat{O}=\hat{O}^{\dagger}$ is Hermitian, then given data $(n_{m},O[n_{m}])$, $m=1...M$ the optimal measurement $\mathcal{O}$ is just the projection of $\hat{O}$ in $\mathcal{S}^{M}_{\mathrm{LL}}=\mathrm{span}\{\rho_{\mathrm{LL}}[n_{1}],\rho_{\mathrm{LL}}[n_{2}]...\rho_{\mathrm{LL}}[n_{M}]\}$. In other words, the reproducing kernel hilbert space (RKHS) of the LLQK [36, 37] consists only of density functionals of the form $O[n]=\mathrm{Tr}\{\rho_{\mathrm{LL}}[n]\hat{O}\}$.

Here we explore the learning abilities of the LLQK in the context of the asymmetric Hubbard dimer. We calculate the LLQK for the dimer using the optimal UCC ansatz states obtained from DC-VQE for specified occupations $n_{1}$ in the relevant interval $\left(0,1\right]$. Since we consider pure state embeddings the fidelities are straightforwardly computed on a statevector simulator using unitary adjoints as $|\langle\Psi(\theta^{\prime})|\Psi(\theta)\rangle|^{2}=|\langle 0|\mathcal{U}^{\dagger}(\theta^{\prime})\mathcal{U}(\theta)|0\rangle|^{2}$. On an actual quantum computer, the kernel entries can be estimated by evolving the initial state $|0\rangle$ with $\mathcal{U}^{\dagger}(\theta^{\prime})\mathcal{U}(\theta)$ and calculating the ratio of all-zero outcomes to the total number of shots. We expect state fidelities to be robust to potential non-uniqueness of the mapping from densities $n_{1}$ to raw UCC parameters $\theta$ as under the Levy-Lieb procedure distances in density space should get mapped to physically meaningful distances in state space. Furthermore, since in this study we fix $2t=1$, we only vary the interaction term $\hat{W}$ by varying $U$ in the dimer Hamiltonian. $U$ can therefore be thought of as a hyperparamter that fixes the kernel and we study kernels $\kappa^{(U)}_{\mathrm{LL}}$ labeled by the $U$ value employed to generate the state embeddings leading to the kernel.

For a first set of learning tasks we consider $F_{\mathrm{LL}}(n_{1})$ as the target function and construct a data set by sampling $n_{1}\in(0,1]$ on a uniform grid with spacing 0.01 and at each $n_{1}$ calculate $F_{\mathrm{LL}}^{(U)}(n_{1})$ and $\Psi^{(U)}(n_{1})$ for different $U$. We then consider two kinds of train-test data splits: (L1) an in-distribution or interpolative learning task and (L2) an out-of-distribution or extrapolative learning task. For L1 we select $M$ training samples distributed as uniformly as possible over the interval $(0,1]$ while including the end points $n_{1}=0^{+}$ and $n_{1}=1$ in the training set. All remaining points are used as test data points. This ensures that the sample mean for $n_{1}$ over the training and test sets is very close and in this 1D example every test data point has at at least one training point on either side. For L2 we select $M$ training samples distributed uniformly over the last $40\%$ of the interval $(0,1]$ and include $n_{1}=1$ in the training set, while all other points are used as test data. In this instance the sample mean of $n_{1}$ over the training and test data is far apart and over 50% of the test data can be reached only through 1D extrapolation. For a given $M$ and specified $U$, we precompute the kernel matrices $\kappa^{(U)}_{\mathrm{LL}}$ for the training and test data sets by evaluating fidelities on a noiseless quantum simulator. For each training set size $M$ considered, we verify the $M\times M$ $\kappa^{(U)}_{\mathrm{LL}}$ matrices are positive semidefinite and then use classical kernel ridge regression (KRR) as implemented in the scikit-learn library [52] to perform fitting and prediction of $F_{\mathrm{LL}}^{(U)}(n_{1})$. The regularization hyperparameter of KRR is set to a small value of $4\times 10^{-15}$ only to ensure numerical stability of the KRR fit process as $M$ approaches the rank of the kernel matrices. Thus we work in a setting where low training fit error is favored. For tasks L1, L2 we also use a classical Gaussian kernel $\kappa_{G}=\mathrm{exp}(-||\mathbf{n}-\mathbf{n}^{\prime}||^{2}/2\sigma^{2})$ to fit and predict $F_{\mathrm{LL}}^{(U)}(n_{1})$ for each U. In order to choose the Gaussian kernel hyperparameter $\sigma$ we use kernal alignment [53, 54, 55] defined for two given kernel matrices $K,K^{\prime}$ as:

Accordingly for each $(U,M)$, we separately optimize $\sigma$ so as to maximize $A(\kappa_{G},\kappa^{(U)}_{\mathrm{LL}})$ on the training data set.

Results for learning tasks L1, L2 as a function of the number of training points $M$ are shown in fig. 4. We use the root mean square (RMS) error on the test data set as our performance metric. Figure 4(a), shows plots for $F_{\mathrm{LL}}(n_{1})$ with $M$=11 training points and corresponding predicted values from task L1 overlaid. The training data points occur at roughly regular intervals on $(0,1]$ and the predictions are indistinguishable from reference data to the naked eye. In fact at $M=11$ the RMS test error for all values of $U$ is under $10^{-3}$ as seen from figure 4(b) where the evolution of the test error is plotted against $M$. We see that with the LLQK the test error quickly falls to under $10^{-2}$ at $M=4$ for all of the $U$ values considered and then continues to fall gradually as more training points are added eventually reaching lower than $10^{-6}$. The RMS test errors obtained by using the Gaussian kernel on task L1 are also plotted in figure 4(b) and show a saturation of the prediction error at around $10^{-2}$ over the range of $M$ investigated. We plot predicted values for $F_{\mathrm{LL}}(n_{1})$ obtained from the Gaussian kernel within task L1 for $M=19$ in figure 4(e). Since the test error is around $10^{-2}$ the fit looks good to the naked eye. In order to understand why the error does not seem to improve beyond $10^{-2}$ we plot the error distribution over the entire test data set on the $n_{1}$ interval $(0,1)$ in figure 4(f). We see that the prediction errors are primarily concentrated near $n_{1}\rightarrow 0$ where $F_{\mathrm{LL}}(n_{1})$ exhibits a sharp increase in its slope [16] and to a lesser extent near $n_{1}\rightarrow 1$. Thus within the chosen setup for task L1, the Gaussian kernel does not generalize as well near $n_{1}\rightarrow 0$ as it is able to elsewhere. Further, on the interpolative learning task L1, we do not find significant improvements in the prediction error of the Gaussian kernel either by increasing the KRR regularization parameter and allowing for higher training loss or by manually tuning $\sigma$. However, since the Gaussian kernel is a universal kernel [56], we expect to be able to reduce the prediction errors by providing more training data near $n_{1}\rightarrow 0$.

In fig 4(c) we plot training and predicted data points as obtained with the LLQK for $F_{\mathrm{LL}}(n_{1})$ within task L2 where the training set is drawn only from the last $40\%$ of the $(0,1]$ interval. We plot the fit obtained for $M=13$ as at this training set size the LLQK RMS test error is under $10^{-3}$ for all U values considered. Figure 4(d) shows the evolution of the LLQK prediction error as a function of $M$ for task L2. We see that in extrapolating $F_{\mathrm{LL}}(n_{1})$ out to densities far away from the training interval the error falls more gradually compared to task L1 but once a sufficient number of training points are provided it declines further to around $10^{-6}$. Thus in this simple model of the asymmetric Hubbard dimer, the LLQK is able to generalize with high accuracy regardless of how the training and test data are distributed over the domain of $n_{1}$. In contrast, the Gaussian kernel is not expected to generalize with naive KRR to extrapolative test data and especially with the small regularization hyperparameter originally chosen, we see from Figure 4(d) that prediction errors are unacceptably high. We find (not shown) that increasing the KRR regularization parameter on task L2 for the Gaussian kernel from $4\times 10^{-15}$ to around $10^{-5}$ improves prediction errors somewhat from $10^{0}$ to around $10^{-1}$ at which level they remain as a function of $M$.

The observed behavior of the LLQK for the dimer suggests it leads to a restricted model with a low effective dimension in feature space. For $M=16$ we show the decay of the singular values of the $\kappa^{(U)}_{\mathrm{LL}}$ matrices for different $U$ in figure 5(a). It is apparent that the magnitudes of the singular values initially decay rapidly suggesting a low effective dimension. Additionally we find that as the size of the training set $M$ increases, the numerically computed rank of the $\kappa^{(U)}_{\mathrm{LL}}$ matrices stops growing around 16 - 21 irrespective of how the training data is distributed over the domain of $n_{1}$. So even as nominally the 4-qubit system is associated with a 256 dimensional feature space, the Levy-Lieb embedding for the dimer is itself restricted to a small subspace and once enough training points are provided to effectively span important features within this subspace, operator measurements fit to the training data are predictive for test data if the function being predicted is of the form $O^{(U)}[n_{m}]=\mathrm{Tr}\{\rho_{\mathrm{LL}}^{(U)}[n_{m}]\hat{O}\}$. To emphasise this latter point we consider the transferability of the $\kappa^{(U)}_{\mathrm{LL}}$ kernels for a specified $U$ in terms of learning density functionals associated with a different interaction strength $U^{\prime}$. To this end, firstly in figure 5(b), we show the relative kernel alignment computed using equation 15, between $\kappa^{(U)}_{\mathrm{LL}}$ matrices associated with different $U$. Not surprisingly the alignment $A(\kappa^{(U^{\prime})}_{\mathrm{LL}},\kappa^{(U)}_{\mathrm{LL}})$ decreases as $|U-U^{\prime}|$ increases. In figure 5(c) we show the prediction errors associated with an interpolative learning task where we attempt to learn the interacting Levy-Lieb functional $F_{\mathrm{LL}}^{(5)}(n_{1})$ using the non-interacting LLQK $\kappa^{(0)}_{\mathrm{LL}}$ and vice versa. We find in this case that the prediction errors level out just above $10^{-2}$ as the number of training data points approaches the embedding rank, which is qualitatively different behavior than what we observed for each kernel $\kappa^{(U)}_{\mathrm{LL}}$ when it was used to learn functionals for the same $U$.

Next, we show that results obtained in terms of learning $F_{\mathrm{LL}}(n_{1})$ are not specific to its particular form and other density functionals are expected to behave similarly. Accordingly we consider two additional functionals the first being the kinetic energy functional

and the second being the so-called Hartree-Exchange-Correlation ($HXC$) functional defined as

The kinetic energy functional $T^{(U)}_{\mathrm{LL}}[n]$ is interesting because even for the dimer, it is a non-monotonic function and its form evolves non-trivially with $U$ (see Fig. 6(a)). The $HXC^{(U)}_{\mathrm{LL}}[n]$ functional is instructive because for the dimer as seen graphically in figure 6(c) it has a relatively simple monotonic behavior as a function of $n_{1}$ but in its definition it involves both the interacting $\Psi_{\mathrm{LL}}^{(U)}[n]$ and non-interacting $\Psi_{\mathrm{LL}}^{(0)}[n]$ states. We consider learning $T^{(U)}_{\mathrm{LL}}(n_{1})$ for the dimer in an extrapolative setting similar to task L2 described previously and find as shown in figure 6(a,b) that in spite of its non-monotonic behavior, $T^{(U)}_{\mathrm{LL}}(n_{1})$ is easily generalized by $\kappa^{(U)}_{\mathrm{LL}}$ to far away test data with high accuracy as the number of training data points approaches the kernel rank. For $HXC^{(U)}_{\mathrm{LL}}(n_{1})$ we consider an interpolative setting similar to task L1 above and attempt to learn $HXC^{(U)}_{\mathrm{LL}}(n_{1})$ using $\kappa^{(U)}_{\mathrm{LL}}$. We also increase the KRR regularization parameter in this instance to $10^{-6}$ to improve generalization behavior. Note that $HXC^{(0)}_{\mathrm{LL}}(n_{1})=0$ for the dimer by definition and is ignored. We find once again that because of the involvement of states with $U=0$, kernels $\kappa^{(U)}_{\mathrm{LL}}$ with $U>0$ show a prediction error profile for $HXC^{(U)}_{\mathrm{LL}}(n_{1})$ that flattens out as a function of $M$ (see Fig. 6(d)) in sharp contrast to the generalization ability apparent for $T^{(U)}_{\mathrm{LL}}(n_{1})$. Furthermore we see that the level at which the error stops improving depends on $U$. This illustrates the selective nature of the RKHS associated with the Levy-Lieb embedding of densities.

Our results for the Hubbard dimer suggest that if somehow we only had access to the Levy-Lieb quantum kernel but not the states themselves, we could still calculate density functionals of observables to high accuracy given sufficient training data. Since the dimer is a small system, assuming we have the correct kernel at hand say in the form of the LLQK, it is easy to provide enough training data to exceed the effective dimension of the embedding and achieve very low prediction error. As we do not conduct system size dependent studies in this work, we do not draw conclusions about the scaling of the effective dimension associated with the Levy-Lieb embedding and thus data requirements for learning density functionals. In general such an analysis would have to be conducted in a context specific manner as based on complexity theoretic arguments [12, 5], we do not expect generic quantum advantage for machine learning the universal functional of DFT. The learning abilities of general quantum models of the form $f(x)=\mathrm{Tr}(\hat{O}\rho(x))$ have been analyzed recently by several authors [5, 35, 36, 57, 37]. Huang et al showed that while quantum kernels of the type $\kappa_{Q}(x,x^{\prime})=|\langle x|x^{\prime}\rangle|^{2}=\mathrm{Tr}(\rho(x)\rho(x^{\prime}))$ can be expected to learn a very general class of functions, in the worst case, the training data requirements to achieve small prediction errors could be very large. Furthermore recent theoretical works [57, 37] have also discussed the intrinsic weakness of fidelity based kernels with regards to exponentially decaying off-diagonal kernel matrix elements in high dimensions and suggested projected kernels based on reduced density matrices (RDMs). Fortunately for most observables relevant to DFT and materials science mapping densities to feature spaces based on RDMs would also suffice [10]. Additionally with regards to classical machine learning, it should be noted that techniques related to machine learning DFT [8, 58, 59, 30, 60] employing both implicit kernel methods and deep learning are now very mature and results from the Gaussian kernel based KRR used in this work for illustration purposes on the asymmetric Hubbard dimer are not meant to suggest quantum advantage for DFT.

We discussed the constrained-search formulation of density functional theory (DFT) within the context of near-term quantum algorithms highlighting the relationship between exact density functionals and underlying wavefunctions. We used parameterized quantum circuits to implement a form of density variational minimization involving run-time calculations of the exact Levy-Lieb functional and illustrated its equivalence to unconstrained wavefunction optimization for obtaining ground states. Interpreting the Levy-Lieb mapping from densities to wavefunctions as a feature space embedding into pure states we demonstrated a quantum kernel that allows us to learn density functionals of observables without obscuring the underlying state space. We hope that our work contributes to improving the explainability of DFT and other reduced variable theories to a broader audience interested in quantum theory.