Learning Density Functionals from Noisy Quantum Data
Emiel Koridon Felix Frohnert Eric Prehn Evert van Nieuwenburg Jordi Tura Stefano Polla




























[e-mail: ]polla@lorentz.leidenuniv.nl





§ ABSTRACT 
    The search for useful applications of noisy intermediate-scale quantum (NISQ) devices in quantum simulation has been hindered by their intrinsic noise and the high costs associated with achieving high accuracy.
    A promising approach to finding utility despite these challenges involves using quantum devices to generate training data for classical machine learning (ML) models.
    In this study, we explore the use of noisy data generated by quantum algorithms in training an ML model to learn a density functional for the Fermi-Hubbard model.
    We benchmark various ML models against exact solutions, demonstrating that a neural-network ML model can successfully generalize from small datasets subject to noise typical of NISQ algorithms.
    The learning procedure can effectively filter out unbiased sampling noise, resulting in a trained model that outperforms any individual training data point.
    Conversely, when trained on data with expressibility and optimization error typical of the variational quantum eigensolver, the model replicates the biases present in the training data.
    The trained models can be applied to solving new problem instances in a Kohn-Sham-like density optimization scheme, benefiting from automatic differentiability and achieving reasonably accurate solutions on most problem instances.
    Our findings suggest a promising pathway for leveraging NISQ devices in practical quantum simulations, highlighting both the potential benefits and the challenges that need to be addressed for successful integration of quantum computing and ML techniques.


Learning Density Functionals from Noisy Quantum Data
    Stefano Polla
    July 2023
====================================================








§ INTRODUCTION


Quantum simulation, NISQ and noise
Following the rapid development of quantum hardware, the last decade saw an explosion in the research on quantum algorithms for noisy intermediate-scale quantum (NISQ) devices <cit.>.
An important goal of this research is identifying avenues towards achieving useful quantum advantage: the application of quantum devices to relevant, classically-intractable problems.
A prominent example is 
the electronic structure problem, which is central to chemistry and material science, and has been a major focus for quantum algorithms due to its inherent quantum nature and significant scientific and commercial importance <cit.>.
Much of the recent research on quantum algorithms for this problem has centered on the variational quantum eigensolver (VQE), a NISQ-tailored algorithm designed to approximate ground states by optimizing a heuristic quantum ansatz <cit.>.
However, several significant challenges hinder the achievement of practical solutions to the electronic structure problem. 
Accurately estimating energies and other properties from a state prepared on a quantum computer is made prohibitively expensive by sampling costs <cit.>.
This is exacerbated by the necessity of applying error mitigation techniques to reduce bias due to circuit noise, which increases the sampling overhead <cit.>. 
Additionally, VQE introduces errors intrinsic to the algorithm, related to the expressibility of the ansatz and the complexity of the optimization landscape.
In practical electronic structure calculations, it is often necessary to solve for multiple system configurations, which increases the overall cost.

Density functional theory and learning functionals
Density Functional Theory (DFT) is the workhorse of modern quantum chemistry and material science, used to investigate the electronic structure of many-body systems at a low computational cost <cit.>.
Central to DFT is the definition of a universal functional, which maps the ground state electronic density to the corresponding kinetic and interaction energy for a given problem family <cit.>.
Although exact universal functionals theoretically exist, their general explicit form is unknown, and their evaluation would necessarily be computationally intractable <cit.>.
Instead, DFT practitioners rely on a vast array of approximate functionals, developed over the decades through mathematical assumptions and physical intuition.
Recently, machine learning (ML) approaches have emerged as powerful tools for designing new DFT functionals  <cit.>.
In particular, deep-learning-based functionals have demonstrated remarkable performance on benchmark problems in chemistry <cit.>.
These models can be trained on mixed datasets generated through various approaches, including DFT based, on other functionals, expensive first-principles methods like coupled-cluster, and experimental data <cit.>.
The synergies between DFT and quantum algorithms have also been explored by a few works, notably in the context of learning functionals on fault-tolerant quantum computers <cit.>, and in using quantum algorithms to supplement results from approximate DFT functionals <cit.>.

classical ML with quantum data
As quantum hardware continues to improve, the prospect of using quantum devices to generate training data that is inaccessible to classical methods is becoming increasingly realistic. 
Such data, derived from the quantum simulation of complex systems, could significantly enhance the capabilities of classical ML models <cit.>.
In turn, these models could generalize from the limited, expensive-to-generate quantum data. 
Moreover, by leveraging the underlying structure, physically-motivated ML models could combine information from all the training data points, filtering out noise and leading to predictions that are more accurate than any individual training point.
While not much research has focused on practically training classical models with quantum-generated data, related concepts exist in the literature. 
In quantum science, ML is often applied to learning properties of a system from measurement outcomes, in tasks like phase classification <cit.>.
Classical shadow tomography <cit.> constructs classical models that can predict many properties of a given quantum state from a limited set of measurements. 
Similarly, recent work demonstrates a technique to extract of classical surrogates from QML models trained on quantum devices  <cit.>. 
Moreover, classical ML models have been employed for noise mitigation in quantum systems <cit.>.

in this paper 
In this work, we explore the application of NISQ devices to generate noisy training data aimed at learning DFT functionals for the Fermi-Hubbard model, with the goal of generalizing and enhancing these datasets.
We focus on limited, noisy datasets of 1,000 data points, which, though small from a machine learning perspective, can already be expensive to generate on a quantum device.
We separately examine the effects of sampling noise typical of NISQ algorithms and the algorithmic errors associated with VQE, analyzing how these factors impact the learned functionals.
The models are benchmarked on both an energy prediction task and a Kohn-Sham-like density optimization task.

The remainder paper is structured as follows: 
Sec. <ref> summarizes the necessary background on DFT, the physical system and the learning task.
Sec. <ref> discusses the dataset generation process, including details on the considered quantum algorithms.
Sec. <ref> describes the machine learning model, training process, and compares various ML regression methods.
Sec. <ref> presents the benchmarking results of the trained ML models.
Finally, Sec. <ref> discusses the findings and provides an outlook for future research. 




§ BACKGROUND




 §.§ Density Functional Theory


Central to DFT are the Hohenberg–Kohn theorems <cit.>, which assert two fundamental principles: (1) the many-body ground state of a system of interacting electrons in an external potential is uniquely determined by its electron density, and (2) the ground state energy can be obtained variationally by minimizing a functional of the electronic density.

To formalize this and set the notation, consider a family of electronic Hamiltonians

    Ĥ = Û + T̂ + V̂

where Û represents a fixed electron-electron interaction, T̂ is a fixed kinetic energy term, and V̂ is a free external potential term. 
The expectation value of the potential ⟨V̂⟩ is defined as a linear functional of the electronic density .

The first Hohenberg-Kohn theorem ensures the existence of the universal functional

    ℱ: _↦ := ⟨ψ_|(Û+T̂)|ψ_⟩,
 
which maps the electron density _ of the ground state |ψ_⟩ of any Ĥ to the respective expectation value of the kinetic and interaction energy .

The functional is termed universal because it is independent of the external potential acting on the system.

The total ground state energy E = ⟨ψ_|Ĥ|ψ_⟩ can be reconstructed from the ground state density as 

    E = ℰ[_] := ℱ[_] + ⟨ψ_|V̂|ψ_⟩,

where ⟨ψ_|V̂|ψ_⟩ is a linear functional of _ by definition.
The second Hohenberg-Kohn theorem states that, for a given V̂, the ground state energy and density of the corresponding Ĥ can be obtained by minimizing this functional: 

    E = min_ℰ[] 
         ,  _ = _ℰ[].


While DFT approaches have been developed to investigate a wide variety of many-body systems, here we focus on a particular lattice model. 
In lattice models, electrons occupy discrete, localized sites, and the electron density  is represented as a vector where each component ρ_j denotes the occupation of the j-th site. 
The functional ℱ[] is the reformulated as a function of this vector, simplifying the computational treatment of the system and making it more amenable to a machine learning scheme. 



 §.§ The Fermi-Hubbard model


In this work, we focus on the Fermi-Hubbard model, a paradigmatic lattice spin-1/2 Fermion model with on-site interaction.
We take this model on a one-dimensional ring of L sites.
The Hamiltonian of this model,


    Ĥ_ = T̂ + Û + V̂_,

comprises a kinetic energy term T̂, ab interaction term Û, and a potential term V̂_. 

All of these can be expressed in terms of fermionic creation and annihilation operators, ĉ_j,σ and ĉ^†_j,σ.
From here on, j ∈{1, ..., L} denotes the site index, and σ={↑, ↓} represents the electron spin.

The number operator is defined as n̂_j, σ = ĉ^†_j,σĉ_j, σ.

The kinetic term describes the hopping of electrons between adjacent sites and is given by:

    T̂ := 
        - t ∑_j^L ∑_σ∈{↑, ↓} (ĉ^†_j,σĉ_j+1, σ + ĉ^†_j+1,σĉ_j, σ),
 
where t is the hopping parameter.
We choose without loss of generality t=1, thereby fixing the units of all energies.
We consider periodic boundary conditions, meaning (ĉ_L+1,σ := ĉ_1, σ).
The interaction term accounts for the on-site repulsion between electrons of opposite spins and is expressed as:


    Û := 
        u ∑_j^L n̂_j, ↑n̂_j, ↓,

where u is the interaction strength. 

The potential V̂_ defines a specific instance of the Hubbard model within the broader family of models characterized by the chain length L and the interaction strength u/t. 
The local chemical potential = (μ_0, ..., μ_L) introduces an additional term at each site to the Hamiltonian:

    V̂_ 
        = ·n 
        := ∑_j^Lμ_j (n̂_j, ↑ + n̂_j, ↓).


Each choice of  defines a unique instance Ĥ_ within the model family, with ground state  | ψ_⟩, and the corresponding ground state energy E.
The elements of the electronic density vector _ are the expected number of particle at each site, summed over the two spin species:

    ρ_j = (n̂_j,↑ + n̂_j,↓)ψ_.


The universal functional ℱ is then a function mapping the vector _ to the scalar = ⟨ψ_|(Û + T̂)|ψ_⟩; we call this the Hubbard functional.
Note that the Hubbard functional depends on the value of the interaction strength u/t, the number of sites L and the geometry of the system.

In our numerical study, we focus on a periodic Hubbard chain of L=8 sites; this system is small enough to allow effective exact diagonalization benchmarks. 
We choose an interaction strength of u/t=4, corresponding to the most challenging regime to simulate when scaling to larger systems <cit.>.
Furthermore, we restrict all solutions to the quarter-filling (L/4 = 2 particles of each spin species) and total-spin singlet subspace, based on the symmetries of the Hubbard model.
A more general functional can in principle be obtained by stitching together functionals defined on different symmetry sectors.




 §.§ The learning task


In recent years, data-driven approaches for learning approximate DFT functionals have emerged, where a machine learning model is trained to map electron densities to energies <cit.>.
For the Hubbard model, Nelson et al. <cit.> demonstrated a method for training a model to learn the functional ℱ using 105,000 density-ground state energy pairs obtained via exact diagonalization. However, the reliance on exact diagonalization limits the feasible system size due to its high computational demands.
In this work, we consider the task of learning an approximation ℱ^ML of the same functional, but from noisy data generated by a quantum algorithm.
Generating such data is expensive, with the cost rising with both the size of the dataset and the desired accuracy.
By teaining a model on a limited dataset, we aim to enable generalization and reduce noise, potentially leading to predictions with greater accuracy than any individual data point.

Our overall approach to learning a DFT functional from noisy quantum data is illustrated in Fig. <ref>: 
The process begins by initializing the Hubbard model with random values of  drawn from a specified distribution (Sec. <ref>).
The ground state problem is then solved for 1000 potentials to obtain a dataset of noisy density-energy pairs (_X, _X), where X indicates the method used to generate the data.
To represent the noise typical of NISQ algorithms, we consider and simulate two methods: expectation value estimation (EVE), which isolates the effect of sampling noise, and the variational quantum eigensolver (VQE), which isolates expressibility and optimization errors.
Exact diagonalization results are used for benchmarking, while the noisy data are employed to train a neural network-based regression model.
This approach allows us to explore how the inherent noise characteristic of quantum data from current near-term devices impacts the training of a machine-learned DFT functional.




§ THE DATASET


Training and testing data sets are constructed by solving configurations of the Hubbard problem with a random on-site potential .
The same configurations are solved with each of the considered methods: exact diagonalization (ED), expectation value estimation (EVE), and variational quantum eigensolver (VQE).
The training set and the test set are generated using 1000 independent random potentials each, the training set will be further split to implement 5-fold cross-validation.
Note that this number of configurations is smaller compared to the 105.000 potentials used in <cit.>. 
We focused on this small-data regime due to the anticipated computational expense of generating data points with quantum algorithms.



Selecting the potentials
To sample the random potential  we follow a modified version of the approach in <cit.>.
For each data point, we first sample a strength parameter W ∈ [0.005t,2.5t] uniformly at random.
We then sample  uniformly at random and calculate its standard deviation
σ() = √(∑_j μ_j^2 - (∑_j μ_j)^2).
If the standard deviation σ() < 0.4 t we accept the potential and add it to our dataset, otherwise we reject and repeat the sampling procedure from the beginning.
This procedure produces a representative distribution of potential of varied strengths avoiding too large energy fluctuations, without requiring the solution of the Hubbard problem instance at this stage.

Exact diagonalization benchmark
We construct an exact dataset to act as a baseline, solving all the problem instances (training and test set) with exact diagonalization.
For every problem instance Ĥ_, we construct the Hamiltonian matrix block corresponding to quarter filling and spin-singlet. 
We diagonalize it, obtaining the ground state |ψ_()⟩ and the energy E(). 
From the state we calculate the expected particle density vector _().
We then compute the kinetic-interaction energy () = E() - ·_().
Thus, from each random configuration , we obtain a pair (_, ) representing the input and output of the exact Hubbard functional ℱ.



 §.§ Expectation value estimation
 

Sampling noise
Assuming the exact ground state |ψ_⟩ of the Fermi-Hubbard Hamiltonian Eq. (<ref>) can be prepared on a quantum device, measuring  and _ will still incur an error called sampling noise.
To reproduce this effect we simulate expectation value estimation, by drawing samples from the measurement outcomes distribution defined by the Born rule on the ground state vector |ψ_⟩ produced by exact diagonalization.

Measurement scheme
To define the measurements for expectation value estimation, we note that all the terms we want to measure are diagonal either in real space (density and Coulomb energy) or in Fourier space (kinetic energy) <cit.>.
To estimate the real-space-diagonal terms, we can sample all the number operators n̂_j, σ at the same time.
Averaging M samples we obtain estimates of each n̂_j,σ, and thus reconstruct an estimate ρ̃_j for the density ρ_j = n̂_j, ↑ + n̂_j, ↓.
The expectation value of the Coulomb energy
⟨Û⟩ = u ∑_j n̂_j,↑n̂_j,↓ 
can be estimated from the same samples by averaging the on-site correlators n̂_j,↑n̂_j,↓.
The kinetic energy operator T̂ is diagonalized by expressing it in terms of the Fourier-space fermionic operators ĉ_k,σ = 1/N ∑_j=0^N-1 e^i j k 2 π/Nĉ_j,σ, 

    T̂ = 2 t ∑_k=0^N-1∑_σ=↑, ↓cos(2 π k/N)  n̂_k,σ

where n̂_k,σ = ĉ^†_k,σĉ_k,σ.
The operators n̂_k,σ can be sampled at the same time; another M samples are taken and averaged to obtain n̂_k,σ, and their results combined to estimate ⟨T̂⟩.
The estimate _EVE of  is finally calculated by summing the kinetic and Coulomb contributions.
For the central limit theorem the estimate of F̃ is unbiased and asymptotically normal, with variance proportional to 1/M.

Implementation on qubits
Assuming we have a qubit-based quantum device and Jordan-Wigner Fermion-to-qubit encoding <cit.>, the real-space samples are simply obtained by measuring all qubits in the computational basis as n̂_j,σ = (1-Z_j,σ)/2 (where Z_j,σ is the Pauli Z-operator on the qubit that encodes the spin-orbital j,σ).
In order to sample the plane wave number operators n̂_k,σ, we need to perform a basis rotation circuit implementing the fast fermionic Fourier transform <cit.>: a circuit of depth O(N) acting on the two spin sectors independently.



 §.§ Variational quantum eigensolver



definiton of VQE
The VQE is a hybrid quantum-classical method that combines a quantum subroutine with a classical optimizer, aiming to estimate the ground state energy of a Hamiltonian  <cit.>. 
The quantum subroutine prepares on a quantum device an ansatz state |ψ(θ)⟩ = U(θ) |ψ_⟩ dependent on a set of classical parameters θ; and measures its expected energy E(θ) = Hψ(θ) through EVE.
This energy is then minimized over the parameters θ using a classical optimizer which calls the quantum subroutine.

VQE errors (expressibility and optimization error)
The VQE is a heuristic method, which means that convergence to the true ground state energy is not guaranteed <cit.>.
The manifold of states ψ(θ) that can be represented by the chosen ansatz is, in general, much smaller than the Hilbert space of the problem.
This is a feature of VQE: restricting the number of parameters is the only way to make the problem feasible for the classical optimizer.
The minimum energy ansatz state will thus only approximate the ground state, with its energy min_θE(θ) guaranteed to be larger than the ground state energy by the variational principle.
The quality of the results depends on the choice of ansatz, the performance of the classical optimizer, and is further affected by hardware noise and sampling noise.

In this work, we focus on isolating the effect of the errors intrinsic to the VQE due to ansatz expressibility and optimization. 
We use a classical state-vector simulator to extract expectation values without sampling noise, which is studied separately through EVE.
Both the energy and the electronic density of the optimal state are affected by the VQE errors.

Ansatz in VQE
Variational ansatz families fall into two main categories: hardware-efficient and physically inspired. Physically inspired ansätze, like the Unitary Coupled Cluster (UCC) <cit.> consider essential system properties but often require high circuit depths, making their implementation challenging for near-term quantum devices. 
Hardware-efficient ansätze instead maximize the number of parameters per circuit layer, but may still require many layers to approximate the ground state and do not take into consideration the problem symmetries.
We will consider two ansätze in this work: the well-known Variational Hamiltonian Ansatz (VHA) <cit.> and the Number Preserving Fabric (NPF) <cit.> ansatz, which combine the advantages of both categories.

VHA
The VHA is one of the earliest proposed ansätze for VQE study of the Fermi-Hubbard model <cit.>, and has found success in hardware experiments <cit.>. 
This ansatz is inspired by the adiabatic algorithm <cit.> and formally equivalent to the quantum alternating operator ansatz (QAOA) <cit.>.
It distinguishes itself from the usual VQE ansätze by incorporating the values of  defining the problem instance into the ansatz defintion.
The circuit structure for the Hubbard model, for parameters θ={θ_i,j}, is

    |ψ{θ_i,j}⟩ = ∏_i=1^p [ e^-i T̂_e θ_i,0 e^-i T̂_o θ_i,1 e^-i V̂_θ_i, 2 e^-i Ûθ_i, 3] |ψ_0⟩,

thus able to represent a first order trotter decomposition in p slices of the adiabatic evolution to H from H^0, commonly taken as the non-interacting Hamiltonian of which |ψ_0⟩ is the ground state. 
In Eq. (<ref>), T̂_e (T̂_o) are all hopping terms on even (odd) sites in Eq. (<ref>).

NPF
The NPF ansatz takes another strategy which is, in contrast to VHA, independent of the problem. 
It composes the two fundamental interactions contained in a spin-preserving two-body Hamiltonian; a spatial orbital rotation and a double excitation. 
The parameterized composition of these operations we call Q(θ, ϕ), and they are arranged in a brick-wall pattern to maximize the gate density, resulting an expressive ansatz with only local gates. For the exact form of the ansatz, see Ref. <cit.>.

Comparison of VHA vs NPF (nr. parameters and optimizer)
While both the VHA and NPF ansatz are physically inspired and preserve key symmetries like total spin and number of particles, we highlight here some differences between the two. The number of parameters in the VHA equals 4 per layer by inspection of Eq. (<ref>).
The NPF ansatz has two parameters per Q-block, of which there are N_q/2-1 per layer. In our case N_q = 2N = 16, resulting 14 parameters per layer. Although the NPF ansatz has more parameters per layer, its higher gate density means it does not necessarily result in deeper circuits.
VQE cost functions are often highly irregular and non-convex, with numerous local minima. Hence, the choice of classical optimizer is crucial <cit.>. The VHA, being a more structured ansatz, may have a more irregular cost landscape, making it prone to getting stuck in local minima. To address this, we use the gradient-free Constrained Optimization by Linear Approximation (COBYLA) optimizer. For the NPF, the abundance of parameters offers more flexibility but less structure, and empirically, following the gradient in this cost landscape tends to work better. Therefore, we employ the Sequential Least Squares Programming (SLSQP) optimizer for the NPF ansatz.



§ LEARNING THE FUNCTIONAL


In this section we present the chosen machine learning models and the details of the training procedure.
Alongside the main neural-network model, for which we provide a detailed result analysis in Sec. <ref>A, we present a comparative study of a set of out-of-the box models in Sec. <ref>.



 §.§ Method
 

CNN model
The main ML model we consider is an adaptation of the convolutional neural network (CNN) model proposed by Nelson et al. <cit.>.
The model consists of a sequence of one 1-dimensional periodic convolutional layer with 8 local features and kernel size 3, two dense layers with 128 features each, and a final dense layer with a single output.

Before being processed by the CNN, the inputs _X are normalized.
Each element ρ_j of the density  is normalized using the mean ρ̅ and standard deviation σ_ρ of the whole training set:

    ρ_j ↦ (ρ_j - ρ̅)/σ_ρ =: input_j,

(where we dropped the method index X for convenience of notation).
Conversely, the output of the model is rescaled by the standard deviation of the training energies σ_ and shifted by their mean :

    output↦σ_·output +  := ℱ[].


Cross-validation
In potential real-world application of our method to noisy data from systems too large for exact benchmarking, the training data becomes the sole source of information for model evaluation. Thus, model selection, hyperparameter optimization and validation must be conducted exclusively using the available (noisy) dataset.
To address this, we use 5-fold cross-validation <cit.>.
This involves partitioning the noisy training dataset into five equally sized subsets (folds). 
During each iteration, one fold is set aside as the validation set, while the remaining four folds are used to train the model.
The performance metrics obtained from each fold are then averaged to provide a comprehensive evaluation of the model performance.


Symmetry augmentation
The target functional ℱ[] is invariant under translations, mirror symmetry, and their combination.
This means that ℱ[] does not change under the transformation

    ρ_j ↦ρ_(± j + k) mod N ,  ∀ k ∈{0, ..., N-1}.

We exploit this to augment the training dataset: from each pair (_X, _X) we construct 16 data points, applying the N=8 shifts and N=8 mirror-and-shifts to the density and copying the energy <cit.>. 
This data augmentation is performed after the cross-validation split, on the training and validation splits separately. 
The data points within the training split are then scrambled before dividing in batches for training.

optimization/Early stopping/keras
The model is trained with the Adam optimizer <cit.>.
To avoid overfitting, we implement an early stopping strategy which monitors the validation loss during training, halting when this loss stops decreasing.
The model and its training are implemented using the TensorFlow-Keras library <cit.>.




 §.§ Alternative model selection





model selection
Given the changes in the underlying statistical structure of the data in this work compared to Ref. <cit.>, we evaluate the suitability of the CNN model by comparing its performance with other machine learning models known for their robust regression capabilities and ability to handle noise and outliers. 
The objective of this section is to determine if any of these alternative models can enhance the learning task with our noisy data.

Our findings indicate that the CNN model is the best performing model.
Nonetheless, we present our comparative study for completeness.

models
The ensemble of models to be evaluated comprises 11 regression models, 
including the convolutional neural network (CNN) model from Ref. <cit.>, three linear regression-based approaches known for their outlier robustness (Huber regression  <cit.>, Theil-Sen regression (TS) <cit.>, and random sample consensus (RANSAC) regression <cit.>), five decision tree models (AdaBoost (AB) <cit.>, random forest (RF) <cit.>, gradient boosting (GB) <cit.>), support vector regression (SVR) <cit.>, nearest neighbor regression (KNN) <cit.>, and XGBoost (XGB) <cit.>.

For the CNN model, we adopted the hyperparameters from Nelson et al. <cit.>, while for the other models, we used standard baseline hyperparameters.

x-validation
To evaluate model performance, we present results for models trained on data generated from both the VQE (using an NPF ansatz with depth d=5) and EVE (M=2000 shots) methods. 
These two instances are chosen to represent the typical characteristics of VQE and EVE datasets, and we focus exclusively on them for clarity and readability.

For benchmarking purposes, we also include performance metrics for models trained on exact data. 
In scenarios where our method is applied to noisy data from systems that are too large for exact benchmarking, the cross-validation of training data will be the only resources available. 
Consequently, model selection and hyperparameter optimization must be based solely on these datasets:

For this, we employ 5-fold cross-validation, as introduced in Sec. <ref>. 

plot overview
The comparative performance of the models is illustrated by the mean and standard deviation of their cross-validation scores, as shown in the top panel of Fig. <ref>.

Additionally, both panels display the baseline mean-squared error (MSE) of the noisy data used for training and validation, indicated by dashed lines.

The bottom panel of Fig. <ref> illustrates the models' performance compared to the exact benchmark on a separate test set, generated with a distinct set of potentials not used during training.
While this information is not directly used in the model selection process, it provides insights into the models' ability to generalize to noiseless data. 

plot interpretation
The cross-validation scores for models trained on EVE data never fall below the baseline error of the training data. 
This is because the validation dataset is also noisy, and with EVE's noise being randomly distributed with a zero mean, the model cannot effectively learn or predict the noise pattern.

Among the models tested, those achieving the lowest validation scores on this dataset, aside from the optimized CNN model, are Support Vector Regression (SVR), Gradient Boosting (GB), Random Forest, and XGBoost (XGB). 
These models also perform the best in predicting noiseless data.

We observe the characteristic performance that a cross-validation score close to the baseline of unbiased noise indicates that the underlying model has learned patterns that extrapolate to some extent to noiseless data. 

The performance of models trained on VQE data differs noticeably when tested against exact data. 
As discussed in Section <ref>, the optimization error inherent in our VQE data can lead to an overestimation of the ground state energy for a given model instance. 
This results in biased noise that affects the training process. 
Consequently, the model inadvertently learns to fit this biased noise, which does not exist in the exact data, leading to a worse performance when testing the generalizability. 

When the MSE between the VQE-generated energies and the exact energies (indicated by the dashed blue line) is less than or equal to the predicted MSE from cross-validation, it suggests that the model has absorbed the bias present in the VQE data. 
Specifically, for models 7-13 trained on VQE data, the validation MSE is smaller than the baseline error of the training data.

This observation implies that (1) there is an (expected) bias in the error of VQE and (2) these models learn to predict this bias.


CNN FTW
Based on the cross validation data, the CNN model from Ref. <cit.> performs best on all datasets.
While gradient boosting performs slightly better on the presented VQE dataset (NPF ansatz with d=5), general performance across VQE datasets is at best comparable to the CNN model.
We thus continue further analysis with the CNN model only.



§ RESULTS





In this section, we benchmark the CNN model trained on noisy datasets, analyzing its performance as a function of both the type and magnitude of noise in the training data.
We test the models on two tasks: energy prediction, discussed in Sec. <ref>, and density optimization, discussed in Sec. <ref>.
All benchmark are conducted using the same test set, consisting of 1,000 new problem instances with random potentials _test.
These are solved with exact diagonalization obtaining a set of exact test densities _test and energies _test.



 §.§ Energy prediction



In the energy prediction task, we assess the model's accuracy by comparing its predictions on the exact densities, ℱ^ML[_test], against the exact energies _test.
Figure <ref> summarizes this comparison for four models, each trained on a datasets constructed with a different method: exact diagonalization, EVE (M=1000), VHA-VQE (d=6), and NPF-VQE (d=6).
The figure also reports the noisy training energies , compared to the respective exact value.

The model trained on exact data serves as a benchmark, representing an upper limit of performance.
Models trained on the same amount of noisy data points are expected to perform less accurately.

The EVE data are spread uniformly around the diagonal line, as sampling noise is unbiased.

The model trained on this data demonstrates improved energy predictions, which cluster closer to the diagonal line compared to the noisy training data.

Conversely, the VQE data suffer from a positive bias due to the limited expressive power of the ansatz.
Variational methods always approximate the ground state energy from above, Ẽ_VQE≥ E. 
As a consequence, F̃_VQE typically (albeit not always) overestimates .
We observe this bias is larger for smaller values of , particularly for the d=6 NPF ansatz.
This effect is likely due to the variation in the strength of , which influences the efficiency of VQE – larger potentials cause stronger localization, thereby increasing kinetic and interaction energy as electrons are more confined. 
Localization also reduces ground state correlations, which can enhance the performance of limited-depth VQEs.

The predictions of the models trained on VQE data are distributed similarly to their respective training data, indicating that the model learns the inherent bias.




MSE
To systematically compare models trained on datasets with a varying noise strengths, we summarize the results of the energy prediction task with the mean squared error (MSE)

    MSE_ℱ^ML = 
        1/test size∑_test set |ℱ^ML[_test] - _test|^2.

This is compared to the raw MSE of the method used to generate the training data, calculated directly from the training set as 

    MSE_X = 1/train size∑_train set |_X - |^2.


EVE-trained models
Figure <ref> shows the performance of the model trained on EVE data as a function of the number of samples M.
The squared error in the training data decreases proportionally to 1/M, consistent with the standard sampling limit.
(The factor 4 arises from the average ground-state variance of the two operators T̂ and Û, which are sampled separately as discussed in Sec. <ref>.)
The model's predictions improve on this MSE by up to factor 20 in the case of large sampling noise.
As the sampling noise decreases, the model's performance saturates to that of the model trained on exact data.
It is important to note that the total number of samples needed to generate the training dataset is 1000 M, uniformly distributed over the 1000 different problem instances.
The improvement is due to the model's ability to synthesize the information from all these training instances, distilling the underlying trend while filtering out the unbiased noise.

VQE-trained models
Figure <ref> presents the same performance benchmark for the model trained on VQE data.
The MSE of the VQE, for both ansätze considered, decreases exponentially with the depth of the ansatz.
Here, the depth d represents the number of repeated layers in the ansatz, while the upper-axis labels indicate the number of parameters.
For depths larger than d=6 the optimization of the VHA becomes significantly harder, thus we explore the larger-parameter-number regime with the less-structured NPF ansatz.
The learned functional generalizes the training dataset successfully, showing a MSE comparable to the respective training data for all depths and both ansätze.
Unlike the models trained on EVE, those trained on VQE data do not manage to learn the underlying functional while filtering out noise. 
This limitation is due to the bias intrinsic to the VQE errors, which is learned by the model.
Although there is some quasi-stochastic error in the VQE results, caused by the optimization converging to different local minima of the energy, the dominant factor appears to be the biased error due to the limited expressibility of the ansatz.





 §.§ Density optimization





Motivation
The main application of DFT functionals lays in solving new problem instances through density optimization.
The goal of density optimization is equivalent to usual Kohn-Sham self-consistent field, but rather than optimizing a set of single-particle orbitals together with an exchange-correlation functional, we directly optimize the total energy functional in Eq. (<ref>) with respect to the density  (subject to appropriate constraints).
In this section we explore the performance in this application of the ML functional learned from noisy data using the CNN model.

Task description
For each test set potential _test and each considered learned functional F^ML, we construct the total energy functional

    ℰ^ML[] = ℱ^ML[] + _test·.

We then minimize ℰ^ML[] under the constraints of ρ_j ∈ [0, 2] for each site j (positive and bounded occupation) and ∑_j ρ_j = 4 (fixed total number of particles); this yields predictions for the ground state density ^*_ and energy E^* = ℰ[^*].
The minimization is performed using the Sequential Least Squares Programming (SLSQP) method implemented in the SciPy package <cit.>, exploiting automatic differentiation of the CNN model to obtain analytical gradients of ℰ[].

Description of Figure <ref> 
Figure <ref> shows the resulting errors, for a selection of models trained on datasets constructed with different methods with different accuracy: exact diagonalization, EVE (M=1000), NPF-VQE (d=6), and NPF-VQE (d=11).
For energy we consider the absolute deviation |E^* - E_test|, and for density we consider the error in ℓ^2-norm 
    ‖^* - _test‖_2
        =
        √(∑_j |ρ^*_j - ρ_test j|^2).

In the scatter plots we show these errors for each of the test set instances versus the strength of the potential  —
as the average value of the potential μ̅= ∑_j μ_j shifts the functional Eq. (<ref>) by a constant and does not affect optimization, the strength of the potential is characterized as its standard deviation sttdev() = ‖ - μ̅‖_2.
To the right of each scatter plot, the distribution of the errors in summarized in a histogram.
As all histograms display bimodal distributions, with two clearly separated lobes at lower- and higher-errors.
On each histogram, we show the line that separates the two modes, and report the fraction of training data in the higher-error mode.

Interpretation of Figure <ref>
The results in the left panel of Fig. <ref> show that, even when using the model trained on exact data, density optimization does not always converge to the correct result.
In particular, for the largest potential strengths the optimization yields densities and and energies with a large error.
We attribute this to the limited availability of training points near the edges of the potential distribution.
In the region of smaller potential strengths, at the center of the distribution of , has more data available; there the model can learn and generalize correctly.
Farther from the center of the distribution, where the potential strength is larger, few data points are available for a wider area of possible cases to cover.
We estimate 16% of the test set potentials are in this edge area; and a similar fraction of training points are available there.
While this unavoidable for a continuous and unbounded distribution of problem instances, one can mitigate the problem by training the model on a set of potentials wider than the application range.

Conversely, for the smallest potentials, the density optimization sometimes yields correct energies and wrong densities.
This occurs because, for very small potentials, the considered Hubbard model can exhibit two quasi-degenerate states with similar energies but significantly different densities.
As a consequence ℰ[] has multiple local minima for very distant densities, and the local optimizer we use sometimes finds a wrong density with a good associated energy.

Both of these effects, for the smallest and largest potential strengths, affect the results of the models trained using noisy data in a qualitatively similar way.
Noisy training data causes an overall increase of error on the densities and energies, as well as shrinks the range of potentials on which the results of energy optimization are reliable.
In the case of VQE with lower depth (NPF d=6), we can see a concentration of the energies due to the large positive bias of the training data, which dominates the energy error.




§ DISCUSSION AND OUTLOOK


conclusion
In this work, we demonstrated the application of machine learning models motivated by density functional theory to generalize and enhance a set of data subject to noise characteristic of near-term quantum algorithms.
We showed that meaningful density functionals can be learned from a small amount of noisy training data and benchmarked the performance of these models in their usual DFT applications.

Comparing the performance of the learned models to the respective training data, we showed a significant improvement is only achieved when the noise in the input data is unbiased (Fig. <ref>).
This implies that learning techniques are well-suited to improve results from quantum algorithms that suffer from sampling noise.
However, even when the trained models learn the inherent dataset bias, the proposed technique could be useful to generalize the dataset to new problem instances (Fig. <ref>).

In the near future, quantum devices might be able to generate approximate solutions for quantum simulation problems that are beyond the reach of classical computation but remain too noisy to be directly useful.
Nonetheless, this data could be valuable for supplementing the training of classical machine learning models, possibly in combination with other classically-generated data.
In this context, it is important to characterize the effect of noise in quantum data on the trained models.

Joint noise study
Our study focused on analyzing the (biased) VQE error and (unbiased) sampling noise separately. 
However, in practical implementations on quantum hardware, both sampling and optimization noise critically influence the performance and accuracy of VQE algorithms. 
Additionally, hardware noise is a significant factor. 
Investigating the combined effects of these noise sources presents a promising direction for future research.

Scaling with number of samples and system size
Furthermore, we concentrated on learning from a very limited dataset, comprising only 1,000 training and validation instances. 
To better understand the trade-offs involved in achieving functionals with useful precision, it would be relevant to study the performance of learning as a function of both dataset size and problem size.

Additionally, advancing towards realistic applications will require the development of techniques that enable learning from mixed datasets, combining results from approximate classical algorithms and quantum data. 
In this context, transfer learning might offer valuable insights  <cit.>.

Which functional to learn
In this study, we focused on learning Hohenberg-Kohn density functionals, a cornerstone of density functional theory with extensive literature support; however, the choice of learning targets can generally vary depending on the specific application.
In the future, we envision betond-classical quantum computations contributing valuable data to further enhance a wider range of existing machine learning-based functionals.

Another potential learning target is the exchange-correlation functional used in Kohn-Sham density functional theory, which differs from our functional by a classical estimate of the kinetic energy expectation value, based on a product state. 
Alternatively, learning targets could be drawn from 1-particle reduced density matrix functional theory, where the underlying compressed representation of the quantum state is its 1-particle reduced density matrix  <cit.>.

Finally, one could consider directly learning the mapping → E, also known as the Hohenberg-Kohn mapping <cit.>. 
A comparative study could help determine which of these learning targets is most resilient to noise in the data, offering better adaptability in real-world scenarios.



§ ACKNOWLEDGEMENTS

We thank Vedran Dunjko and Patrick Emonts for their useful feedback.
This work was supported by the Dutch National Growth Fund (NGF), as part of the Quantum Delta NL programme.
S.P. and E.K. acknowledge support from Shell Global Solutions BV.
J.T. acknowledges the support received from
the European Union’s Horizon Europe research and innovation programme through the ERC StG FINE-TEA-SQUAD (Grant No. 101040729).
This publication is part of the `Quantum Inspire - the Dutch Quantum Computer in the Cloud' project (with project number NWA.1292.19.194]) of the
NWA research program `Research on Routes by Consortia (ORC)', which is funded by the Netherlands Organization for Scientific Research (NWO).
The views and opinions expressed here are solely those
of the authors and do not necessarily reflect those of the
funding institutions.
Neither of the funding institutions
can be held responsible for them.



§ CODE AND DATA AVAILABILITY

The density-energy pair data used for training the machine learning models, along with the code used to generate it and to reproduce the plots in this paper, is available at the repository: https://github.com/StefanoPolla/DFTQMLhttps://github.com/StefanoPolla/DFTQML.



§ COMPETING INTERESTS

The authors declare no competing interests.



§ CONTRIBUTIONS

E.P., E.K. and S.P. conceived the project and developed a first prototype of the code.
S.P. developed the codebase for data generation and learning.
E.K. developed and optimized VQE.
F.F. performed the model selection studies.
E.K., F.F. and S.P. wrote the manuscript, which all authors contributed to reviewing.
J.T., E.vN. and S.P. supervised the work and provided funding.



