{'fraction_non_alphanumeric': 0.0321866485013624, 'fraction_numerical': 0.011126248864668483, 'mean_word_length': 4.68513339070568, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reliable emulation of complex functionals by active learning with error control

A statistical emulator can be used as a surrogate of complex physics-based calculations to drastically reduce the computational cost. Its successful implementation hinges on an accurate representation of the nonlinear response surface with a high-dimensional input space. Conventional ``space-filling'' designs, including random sampling and Latin hypercube sampling, become inefficient as the dimensionality of the input variables increases, and the predictive accuracy of the emulator can degrade substantially for a test input distant from the training input set. To address this fundamental challenge, we develop a reliable emulator for predicting complex functionals by active learning with error control (ALEC). The algorithm is applicable to infinite-dimensional mapping with high-fidelity predictions and a controlled predictive error. The computational efficiency has been demonstrated by emulating the classical density functional theory (cDFT) calculations, a statistical-mechanical method widely used in modeling the equilibrium properties of complex molecular systems. We show that ALEC is much more accurate than conventional emulators based on the Gaussian processes with ``space-filling'' designs and alternative active learning methods. Besides, it is computationally more efficient than direct cDFT calculations. ALEC can be a reliable building block for emulating expensive functionals owing to its minimal computational cost, controllable predictive error, and fully automatic features.


I. INTRODUCTION
Theoretical research in chemistry and materials science is often hampered by computational bottlenecks in applying complex physical models to predict the microscopic structure and physicochemical properties of many-body systems. Statistical and machine learning (ML) methods, such as Gaussian process (GP) regression, basis expansion methods, and neural network models, have been widely used as a surrogate of complex physical models to emulate the atomic force fields, density functionals, chemical reactions, and diverse properties of materials and chemical systems [1][2][3][4][5][6][7][8] . Here a typical aim of a ML approach is to learn the map with a high dimensional input or a functional input from observations. Trained on a set of pre-specified inputs, a statistical emulator can obtain accurate predictions when a test point is close to the training inputs [9][10][11] .
Constructing a statistical emulator often starts with selecting a set of "space-filling" samples in the input space, such as uniform samples or Latin hypercube samples (LHS) 12 . The statistical emulator will be trained based on the outputs of the simulation at these pre-selected inputs. The idea is to ensure that for each test point in the input space, there are a sufficiently large number of inputs near this test point, where the outcomes (such as the potential energy or atomic force in emulating the first-principles calculations) were observed. Using the space-filling samples and assuming a set of regularity conditions, the predictive error of some ML methods, such as the GP regression, is guaranteed to converge to zero, when the sample size goes to infinity 13 . The spacefilling samples have two main drawbacks in emulating physics-based calculations such as molecular dynamics simulation (MD). First, the input values, such as atomic positions or pairwise inter-atomic distances in MD, may not be directly controlled. Second, physics-based modeling of molecules and materials often involves a high or infinite dimensional input space, such as the external potential function and atomic configurations. If the statistical emulator is trained in a small parameter subspace, the predictive accuracy of the emulator can be dramatically degraded outside the subspace where no training input is available. This scenario often occurs when the outcome is required in a slightly different system or at another experimental condition, such as different temperatures or pressures, whereas the statistical emulator becomes inaccurate if it is only trained on one set of conditions. Obtaining accurate predictions of any test input with a controlled predictive error is important to bypass expensive computations based on complex physical models.
To address the fundamental difficulty in predicting functions and functionals with a high-dimensional input space, active learning has become one of the most promising techniques to sequentially reinforce the emulation of physics-based modeling 14 . The active learning approach can be implemented in two broad scenarios. The first one is an online scenario or "on-the-fly" prediction, where the test inputs sequentially come and one does not know what future input needs to be predicted [15][16][17] . The online prediction scenario has many applications, such as Bayesian optimization 18,19 and model calibration 20,21 , where the outcome of the simulation needs to be sequentially predicted for every sampled input. The second scenario is to sequentially select the samples to run a computer simulation for predicting the entire input space, where the uncertainty of the emulator is often used for selecting the next design run 22 . In this work, we focus on the online prediction scenario and test our approach using the classical density functional theory (cDFT) calculations for one-dimensional (1D) hard-rod systems, with different external potentials. The availability of exact densities and free-energy functional for various 1D systems of hard rods allows us to validate the efficiency of the surrogate model with the ground truth. Here the fundamental difficulty is on a user-specified external potential-an arbitrary functional input with infinite dimensions. For constructing statistical emulators, we assume that no parametric form of this functional input is available. We demonstrate that the new integrated approach provides a high-fidelity prediction of particle density profiles and grand potentials with a controlled error in the probabilistic sense, more accurate than a GP emulator with conventional "space-filling" designs, and it is more computationally scalable than direct cDFT calculations.
We highlight the key novelty of this work in developing the surrogate model and error control criteria for functional mapping. First, while some popular criteria, such as the D-optimality, were introduced in the active learning framework 15,23 , the threshold of determining whether a test input can be reliably predicted may be hard to choose, as the threshold may not be directly related to the error in prediction. In practice, one may need to tune the threshold for different systems for a given error bound. Here, we use the internal uncertainty assessment of the GP emulator to decide whether a test input can be predicted precisely. For any given threshold of predictive error and probability tolerance bound, our approach automatically defines the criterion to control the predictive error below the threshold with a large probability satisfying the probability tolerance bound (see Section III B). We demonstrate that our approach can identify whether an upcoming input can be reliably predicted, and thus it has much better performance than the conventional statistical emulation by training the emulator with a spacefilling design, such as random sampling. Furthermore, we derive an iterative formula that reduces the computational complexity for Gaussian process models with augmented samples. In each iteration, updating our algorithm only requires O(n 2 ) operations, while a direct approach requires O(n 3 ) operations. Other emulation approaches, such as those enforcing the invariance symmetries, can be easily included in our integrated framework for reliable predictions with a controlled error.
The remainder of the paper is structured as follows: In Section II, we provide a brief overview of the clas-sical density functional theory in the context of the 1dimensional (1D) hard-rod model used in this work. In Section III A, we introduce the parallel partial Gaussian process for emulating physics-based calculations with vectorized output on massive grid points, such as particle densities from cDFT calculations. Then, in Section III B, we introduce our proposed method to control the predictive error. The computationally scalable update of the algorithm and its computational cost are discussed in Sections III C and III D, respectively. The numerical comparison in Section IV demonstrates that our approach outperforms the conventional GP emulators with random-sample designs and the active learning approach with the D-optimality criterion. We conclude the paper along with perspectives on some future research directions in Section V.

II. CLASSICAL DENSITY FUNCTIONAL THEORY
Classical density functional theory (cDFT) is a statistical-mechanical method that describes the thermodynamic properties of equilibrium systems in terms of the density profiles of individual particles 24,25 . In a nutshell, cDFT calculations are based on the minimization of a grand potential functional that, for systems consisting of only one type of particle, can be expressed as where ρ(r) denotes the particle density at position r, F id [ρ] is the free-energy functional of an ideal gas, and F ex [ρ] is called the excess free-energy functional, V ext (r) denotes a one-body external potential, and µ is the chemical potential of the particles. The ideal-gas functional has an exact form where Λ is the thermal wavelength and β = 1/(k B T ) with k B being the Boltzmann constant and T the absolute temperature. While T and µ are thermodynamic variables that define the equilibrium condition, Λ and F ex [ρ] depend on the intrinsic properties of the system such as the particle mass and inter-particle potential. In this paper, we set β = Λ = 1 as the units of energy and length.
One key premise of cDFT calculations is that V ext (r) is uniquely determined by ρ(r) such that, given an analytical expression for F ex [ρ], ρ(r) can be solved by minimizing Ω[ρ] and, subsequently, all thermodynamic properties of the system can be derived. The equilibrium density profile minimizes Ω and thus satisfies the Euler-Lagrange equation: Eq.
(3) provides a functional map between the one-body density ρ(r) and the one-body external potential V ext (r). Because in general F ex is a complicated functional of ρ(r), solving ρ(r) from the Euler-Lagrange equation is often computationally demanding. Besides, the formulation of an accurate expression for F ex is a formidable task for complex molecular systems. In this work, we demonstrate that statistical emulation will be helpful to solve both problems.
For proof of concept, we consider a statistical emulator of cDFT for one-dimensional (1D) systems consisting of identical hard rods (HR). For such systems, the excess free-energy functional, thus the grand potential, is exactly known 26 . The analytical results were derived decades ago and are reproduced in Appendix A. Given any external potential V ext (r) and chemical potential µ, the density profile of hard rods can be numerically solved (e.g., by Picard iteration) from Eq. (3), and it is plugged into Eqs. (3), (A2), and (A3) for computing the grand potential and other thermodynamic quantities. For a given system defined by µ and V ext (r), statistical emulation aims to predict the particle density and free-energy functionals as in direct cDFT calculations (illustrated in Fig. 1(a)). Previous work in cDFT 5,27,28 and Kohn-Sham density functional theory (KS-DFT) 29,30 demonstrate the performance of ML models with a single parametric form for the external potential with different choices of parameters. In this work, we show that the active learning procedure is applicable to multiple classes of external potential similar to cDFT calculations with the same intrinsic Helmholtz energy functional. We demonstrate that, when presented with a new class of external potentials, the model can discern whether this new input can be accurately predicted. As mentioned above, functional mapping represents one of the biggest obstacles in statistical emulation as well as machine learning. Filling the entire functional input space directly with a limited number of simulation runs, such as random sampling or the LHS, is not generally possible. When the external potential function is away from the training input set, both statistical emulators and machine learning methods become inaccurate. In this work, we provide a so-lution to this fundamental problem by integrating direct cDFT calculations with a Gaussian process (GP) emulator through an active-learning framework. We derive a probabilistic bound to control the overall prediction error of the particle density based on the internal assessment of the uncertainty from the statistical emulator, making our approach applicable to functional learning in infinite dimensional spaces, as well as reducing the computational costs for inputs that can be accurately predicted. Here our algorithm can reliably learn multiple classes of external potentials, while earlier applications of ML methods in both KS-DFT (e.g. 30 ) and cDFT (e.g. 5,28 ) only demonstrate the performance of their models on one particular class of the external potential. Besides, the active learning scheme developed in this work enables the ML approaches to discern whether outcomes of a test input can be reliably predicted. This ability allows us to control predictive errors as well as substantially reduce the computational cost.
In Fig. 1(b) and (c), the input µ − V ext (.) and output ρ(.) of the statistical emulator are projected onto two t-distributed Stochastic Neighborhood Embedding (t-SNE) factors as commonly used for visualizing highdimensional data 31 . It is evident that the inputs generated from different groups of input functions form different clusters. Similar patterns can be observed for the corresponding particle densities. In addition, group 1 forms a curve in Fig. 1(b) due to the fact that the only change of the input in group 1 is µ, which only results in a shift in the input space, making group 1 mapping the simplest to learn for the statistical emulator, as to be demonstrated later in the section on numerical results. We noticed several points in groups 2 and 3 are close to those in group 1. These points correspond to samples whose input parameters are close to 0 and are thus nearly equivalent to group 1's inputs. These patterns demonstrate the input-output relationship can be captured by distance metrics, as the smaller distance between chemical potentials leads to the higher similarity between the particle densities. This validates the use of kernels in  (4). The different parameters are specified as γ = 0.2, σ = 2, µ = 0 (black curves); γ = 0.25, σ = 1.5, µ = 3 (red curves); γ = 0.33, σ = 1, µ = −2 (green curves); γ = 1, σ = 0.5, µ = 0 (blue curves). (b) The truth (black curve), the predictive mean (blue curve), the 95% predictive credible interval (grey region), and 10 posterior predictive samples (colored dashed lines) after observing 12 data points (black circles).
the Gaussian process to map the distance between inputs (chemical and external potentials) to the correlation between output (particle densities), as the smaller distance in the input space indicates a larger correlation or similarity between outputs. Unlike the dimension reduction approach shown in Fig. 1, we will model the entire density profile jointly in a statistical emulator to encode all information.

A. Parallel partial Gaussian process emulation for vectorized outputs
For demonstration purposes, our statistical emulator is built upon a parallel partial Gaussian Process (PP-GP) emulation approach 32 for representing the map between the particle density, a vectorized output on massive grid points, and the input function from the chemical and external potentials. Other emulators that produce uncertainty quantification for predictions can be used as well. The GP emulator is one of the most commonly used surrogate models for expensive physics-based calculations, and it has been widely used for emulating molecular sim-ulation, the potential energy surface, and atomic force fields 1,7,[33][34][35] . In combination with the dimension reduction technique, the GP emulator has also been used as a surrogate model for approximating the KS-DFT calculation of electron density 30 . Fig. 2 illustrates the Bayesian updating formulation in GP from prior to posterior with one-dimensional inputs and outputs. As shown in Fig. 2(a), the GP prior is a very flexible class of functions, that include smooth or wiggly functions with large or small fluctuation, using different choices of parameters. After assuming 12 Latin hypercube samples from a function f (x) = sin(2πx/10) + 0.2 sin(2πx/2.5), in Fig. 2(b), we plot the predictive mean from GP (blue curve) and the actual value of f (x) (black curve), as well as 10 predictive posterior samples that fit the data points with optimized parameters and the 95% predictive interval (gray area), using a GP emulator implemented in 36 , where the parameters are estimated by the posterior mode. First, we found a small number of samples can constrain the predictive mean and posterior samples to be around the truth. Second, the truth are almost all covered by the 95% predictive interval, indicating appropriate quantification of the predictive uncertainty. This internal uncertainty assessment by the GP emulator enables us to develop the criterion to control the predictive error for emulating functions with high-dimensional inputs in a probabilistic way, which will be illustrated in section III B.
We use the PP-GP emulator to model the particle density in cDFT. To highlight the salient feature of the PP-GP emulator, let x = µ − V ext (.) denote the functional input discretized at p spatial grid points. For any input x, the particle density is defined at k grid points, written as ρ(x) = [ρ 1 (x), ρ 2 (x), ..., ρ k (x)]. The number of spatial grid points used in discretizing external potential and density can be different but, for simplicity, here we let p = k. Given any n cDFT calculations with inputs {x 1 , ..., x n }, the density vector at any spatial grid point is assumed to follow a multivariate normal distribu- is a variance parameter, and R is an n × n correlation matrix between density at any grid point over n input scenarios, with (i, j) entry parameterized by a kernel func- At any input x, the mean of the output is commonly is a row vector of mean basis function. Here we only use the intercept, i.e. h(x) = 1 and q = 1. Note that the mean parameter and variance parameter are distinct for each spatial grid point of the particle density. Having different mean and variance parameters makes the model more flexible to capture the variability of particle densities at different grid points.
We assume that the kernel function is isotropic, i.e., it is a function of the Euclidean distance between any pair of inputs. The correlation is commonly modeled by the power exponential kernel functions or the Matérn kernel functions 37 . In this work, we use the Matérn kernel function with roughness parameter 5/2 to model the correlation between input x a and x b : where d = ||x a − x b ||, with || · || denoting the Euclidean distance, and γ is a range parameter that can be estimated by the training data. The sample path of the GP with the Matérn kernel function in (4) is twice mean differentiable 37 , assuring good predictive performance in emulating smooth functions. Assume we have run cDFT calculations at n different input scenarios {x 1 , ..., x n }, leading to a k × n particle density matrix ρ = [ρ 1 , ..., ρ n ]. We compute the predictive distribution after integrating out the mean and variance parameters for making predictions with the reference prior of the mean and variance parameters 38 , π(θ j , σ 2 j ) ∝ 1/σ 2 j , for spatial grid point j = 1, ..., k. Conditional on estimated range parameterγ and output density at the jth coordinate ρ j , the predictive distribution of ρ j (x * ) at any new x * and coordinate j, follows a noncentered scaled Student's t-distribution with n−q degrees of freedom: wherê with Here the n × q input mean basis matrix fol- can be found in the supplement materials. In PP-GP, the distribution of the output is assumed to be independent at two grid points, and thus we have p(ρ j (x * ) | γ, ρ) = p(ρ j (x * ) |γ, ρ j ), for any j and x * , which allows us to use Eq. (5) for computing the posterior predictive distribution. For further justification of the assumption, see Theorem 6.1 in 32 . The PP-GP emulator was implemented in the "RobustGaSP" package 36 , where the marginal posterior mode estimator is used as a default method to estimate the range parameter γ 32 .
In comparison with alternative methods, the PP-GP emulator discussed above has advantages in terms of computational scalability and internal uncertainty assessment. First, constructing the covariance matrix requires O(n 2 p) operations, while computing the predictive mean in Eq. (6) only takes O(nk)+O(n 3 ) operations, where n, p, and k denote the number of training simulation runs, the number of discretization points in the input functions and particle densities, respectively. The O(n 3 ) operations come from the Cholesky decomposition for computing the inversion and determinant of the covariance matrix. In contrast, building a separate emulator independently of particle densities at each grid point takes O(n 3 k) operations for computing the predictive mean, due to the inversion of the distinct covariance matrices at each spatial grid point. Approaches that model the spatial covariance matrix generally take O(k 3 ) operations, which is even slower as the number of spatial grid points k of the density is large. Second, the uncertainty in predictions can be directly assessed as any quantiles and percentiles of the predictive distributions in (5) have closed-form expressions. The assessed uncertainty allows us to control the predictive error through integrating a numerical solver and an emulator into the algorithm, discussed below.

B. ALEC: active learning with error control approach
Active learning is a sequential design approach in statistics. In an "on-the-fly" online prediction scenario, there is no information regarding the future input to be predicted. If a configuration is new to the current training set, i.e., it is located in a sparsely sampled region, we may not be able to accurately predict this sample by the emulator. For a large functional input space, it is almost inevitable to extrapolate the sampled input space, as the number of simulated runs that populate this input space is sparse. Therefore, it is necessary to establish a criterion for determining whether an upcoming configuration can be reliably predicted.
In this work, we develop the active learning with error control (ALEC) emulator that can automatically control predictive error for simulations with high or infinite dimensional input space. A few criteria were studied before to detect the test case hard to be predicted, such as the D-optimality criterion 15 and the predictive variance criterion 17 . Yet, the threshold of these strategies often requires to be chosen empirically, as they may not be di-rectly related to the predictive error. In our automatic approach, a new criterion has been established such that the threshold can be chosen automatically for any given predictive error bound and probability tolerance bound, based on the internal uncertainty assessment of the emulator. Fig. 3 shows an overview of our strategy. For a given testing input x * , if the criterion that controls the predictive error is satisfied as shown in the yellow box, the particle density profile will be predicted; otherwise, a numerical solver will be called for computing the particle density corresponding to x * , and this new information will be added to the training set and used to update the PP-GP emulator.
Let us first consider controlling the predictive error for particle density at a fixed coordinate j based on the PP-GP algorithm. Given a prespecified error bound δ j > 0, we aim to have the predictive error of a test input x * smaller than δ j > 0 for more than 1−α of the predictions, where α is a small value (e.g. α = 0.05). This means i.e., the probability that the absolute error between the predicted and true densities at the jth coordinate exceeds δ j is less than α. Since small predictive variance from the emulator indicates high-fidelity in prediction, we predict the test input x * if where t α/2 (n−q) is the upper α/2 quantile of a Student's t-distribution with n − q degrees of freedom, K * * andσ 2 j are given in Eqs. (7) and (9), respectively. Otherwise, we will compute this particle density directly from cDFT using a numerical solver, as the assessed uncertainty for predicting this test input is large. This strategy satisfies the probabilistic error control outlined in Eq. (10), as shown in Appendix C.
We have a few remarks regarding this strategy. First, the error bound δ j and probability tolerance α have interpretations. Suppose, for instance, one expects to see the predictive error in more than 95% of predictions smaller than 0.01, one can then let δ = 0.01 and α = 0.05. Second, in the active learning strategy, one reduces the sample space from the original space X to X F,j = {x * : σ 2 j K * * ≤ δ j /t α/2 (n − q)}, where the subscript F denotes the input set that can be reliably predicted. A numerical solver for cDFT needs to be called for any input x * that does not satisfy Eq. (11). Though selecting a smaller δ or α leads to a smaller predictive error, more evaluations from the numerical solver would be required, which increases the computational cost. Thus, a balance between the size of the cut-off value and computational time is needed in practice. Third, the degree of freedom of the Student's t-distribution increases with the sample size. When the sample size is sufficiently large, the t-distribution can be accurately approximated by a standard normal distribution. For instance, when n−q = 300 and α = 0.05, t α/2 (n−q) ≈ 1.968 and the upper α/2 quantile of a normal distribution Z α/2 ≈ 1.960, which only differs from that of a t-distribution by around 0.5%. Thus one can use the normal approximation if the sample size is large. Finally, the error control from Eq. (11) means that if the statistical emulator is a correct representation of the reality, and the number of test samples is infinitely large, we have at least 1 − α probability such that the absolute predictive error is less than δ. This does not preclude the cases where in slightly more than 1−α of the predictions, the predictive error is larger than δ, due to model misspecification or the limited number of test samples. Fortunately, the predictive error at these inputs is typically close to δ, as the assessed uncertainty of these inputs is smaller than the threshold.
We can extend the strategy for controlling the error of the density at a fixed coordinate j to the overall error of density prediction. Consider that we use a uniform cutoff value δ, and that we make predictions for a test input x * if each coordinate satisfies Eq. (11), i.e., We call the criterion in Eq. (12) the maximum variance selection criterion. Under this criterion, we are able to derive the probabilistic bound to control the predictive error of particle density at each coordinate, and the root of mean squared error (RMSE), defined as The detailed derivation of Eq. (13) and (14) is given in Appendix C. Eqs. (13) and (14) imply that the overall predictive error can be controlled in a probabilistic way under the maximum variance selection criterion in Eq. (12). In practice, the cut-off in (12) is too conservative, and we rarely see more than α samples larger than the claimed threshold. A less restrictive criterion to control the predictive error is by predicting any input set x * such that We call the criterion in (15) the average variance selection criterion. Eq. (15) controls the predictive error of vectorized outputs (particle densities) in the average sense, while it does not guarantee the error control as in Eq. (12). In practice, one may select a large number of points for numerical solvers to run using criterion (12), thereby requiring a larger computational cost, while the predictive error control can be achieved using criterion (15) with much less training samples. To quantify whether the error is indeed controlled based on our threshold, we plot the empirical cumulative distribution function (CDF) of the absolute errors in Fig. 4 using the criterion (15) with different error bounds and probability tolerance thresholds α for all 4 classes of input functions considered herein. In the upper middle panel, for instance, the vertical solid line marks the chosen threshold δ = 0.01 and the purple dashed line labels the percentile corresponding to 1 − α = 0.95. On the x-axis, we mark the values of the 95th percentile of the absolute errors for each group of input, all of which are less than 0.01. In the remaining panels of Fig. 4, the empirical CDFs of the absolute errors for various α and δ values are illustrated. In general, the smaller error bound or probability tolerance threshold, the smaller the predictive error is, while more samples are required for achieving this accuracy level. Thus, a balance between the computational cost and accuracy level is generally required.

C. Sequential update of the ALEC emulator
When calculating the predictive mean and predictive correlation in Eqs. (6)-(7), one needs to solve R −1 z, for an n-dimensional real-valued vector z ∈ R n . In the PP-GP emulator 32 , we implement the Cholesky decomposition of the correlation matrix R = LL T , and use the forward and backward substitution algorithms to compute R −1 z. In the ALEC emulator, whenever a sample x * ∈ R p is added to the previous training set with n samples, the PP-GP model needs to be updated accordingly. Directly applying this approach will require O(n 3 ) oper-ations for computing the Cholesky decomposition each time when the sample size is n. To reduce the computational cost, we introduce a new approach that takes only O(n 2 ) operations in performing the Cholesky decomposition.
To update L with the addition of x * , we denote the correlation matrix of previous n design elements as R n = L n L T n , with L n being the Cholesky decomposition of R n . Next, we write the updated correlation matrix R n+1 : where r n (x * ) = [K(x * , x 1 ), ..., K(x * , x n )] T . In evaluating R n+1 = L n+1 L T n+1 , it is clear that the first n rows of the lower triangular part of L n+1 is identical to the lower triangular part of L n L n+1 = L n 0 n v n l n .
Therefore, only the last row of L n+1 needs to be computed: with v n (j) and r n (j) denoting the jth element of v n and r n (x * ), respectively, and L n (i, j) denoting the (i, j) entry of L n . Eqs. (16) and (17) indicate that updating L n+1 requires O(n 2 ) computational operations, for i, j = 1, ..., n.
After obtaining the Cholesky decomposition L n+1 we can update other terms required in computing the predictive mean and variance.

D. Algorithm and the computational cost
We summarize the ALEC emulator in Algorithm 1. Given the initial PP-GP emulator, predictive error bound δ, and probability tolerance threshold α, the ALEC emulator determines whether the particle densities can be reliably predicted for a new test input. The output contains the predictive particle densities for all test inputs and they are plugged into the energy functionals to compute the energy. Note that the range parameter in the kernel needs to be numerically estimated, which has the largest computational cost in the algorithm. When more training samples are added, the previously estimated range parameterγ may no longer be accurate if our initial model has a small training sample size. To compromise the computational cost and accuracy in predictions, the range parameter in the ALEC emulator is re-estimated for every 50 training samples, if the training size is less than 350. When the training size is large, the estimated range parameter does not change much. In general, the frequency of re-estimating the kernel parameter depends on the computational budget.
Denote the initial training size as n ini and the number of augmented samples as n aug . The computational cost of constructing the Cholesky decomposition of the initial ALEC emulator is O(n 3 ini ), and the total computational order for updating the Cholesky decomposition is then O( The computational complexity of the online updating approach reduces the cost compared to directing computing the Cholesky decomposition at each iteration, which costs O( nini+naug n=nini+1 n 3 ). The computational complexity is significantly faster than the time required to calculate the density using a numerical solver for computing the particle density at each test input. Indeed, the largest computational cost of our approach comes from the numerical solver for the training input set, as will be seen in Fig. 8. Since only a small fraction of samples are used as the initial training and augmented sample, the computational cost of the ALEC emulator is much smaller than running numerical solvers for each density.

8:
Using the numerical solver to compute the corresponding output ρ(x * ); 9: if the training size is a multiple of 50 and the training size ≤ 350 then end if 16: end for Output: The number of augmented size naug, predicted mean, and variance of particle densities for remaining nt − naug samples.

IV. NUMERICAL RESULTS
Here we evaluate the predictive performance of our surrogate model based on particle densities generated from four classes of external potentials (closed-forms shown in Appendix B) with the length of hard rods a = 1 and the domain size L = 9. The values of chemical external input and parameters in each external potential are changed when generating the data. We assume the external potential functions and chemical potentials are used as inputs, while the parametric forms of the external potential functions are not known. We will test the predictive performance separately for each class of external potential functions, and jointly by combining all classes of external potential functions, in Section IV A and Section IV B, respectively. Section IV C gives the predictive performance of our approach for predicting a new class of densities.
We record the out-of-sample root mean squared error (RMSE) of particle density and the grand potential of the test samples from is the predicted density of ith hold out sample at jth coordinate, and Ω(x * i ) andΩ(x * i ) are the predicted and true grand potential for ith hold out sample, respectively. Other criteria, such as the proportion of the test samples covered in the 95% predictive intervals and the average length of the predictive intervals, are provided in supplementary materials. Note that we only compute the predictive error on those inputs we predict, not including those initial samples or the augmented samples computed by the numerical solver.
We generate 2000 samples for each group with LHS, and compare our surrogate model with three other commonly used approaches. The first two approaches use conventional random sample (RS) designs, where we randomly draw n train samples from the LHS samples as the training dataset and use the remaining samples as testing data. We have two distinct training sample sizes for RS samples. The sample size in the first RS sample (RS1) is specified as the initial training size in ALEC, while the sample size in the second RS sample (RS2) is specified as the final training size in ALEC. For predicting any test input, the training sample size in active learning is neither smaller than the training sample size in RS1 nor larger than the training sample size in RS1. We have also implemented the active learning algorithm with Doptimality 15 , where the details are given in Appendix D.
A. Particle density and grand potential prediction for each group of external potential field separately We first test different approaches for four families of functional inputs separately. The density profiles are initially trained using PP-GP with n ini = 20 inputs per group. Then following our ALEC algorithm, the surrogate model is sequentially updated by augmenting the train data with samples selected by the average variance selection criterion in Eq. (15). After obtaining the predicted densityρ for the remaining test samples, we predict the grand potential Ω by pluggingρ into the energy functional in Eq. (A2) in Appendix A. Figure 5 displays the overall performance of different surrogate models (more detailed results can be found in the table in the supplementary materials). Here the threshold in the ALEC algorithm is set as δ = 0.01 and α = 0.05, representing a strategy for having less than 5% of the coordinate to have an absolute error of density prediction larger than 0.01. The first group of external potentials has the least flexibility, as the only change in the input is µ. Thus, predicting the density profile for FIG. 5. The bar plots and line plots show the out-of-sample RMSE of density ρ, and RMSE of grand potential Ω, respectively. The ALEC emulator is compared with other methods with three other designs: RS1, RS2, and D-opt (Active learning with D-optimality). For the ALEC emulator, we use δ = 0.01 and α = 0.05 as the threshold, and the augmented sizes for groups 1-4 are 0, 7, 21, and 504, respectively. In Dopt, the augmented sizes for groups 1-4 are 0, 7, 21, and 582, respectively. The number of training sample sizes from RS1 and RS2 are the same as the initial sample size and terminal sample size in ALEC.
the first group of functions is the easiest. Models with only 20 inputs in the conventional RS design can predict particle densities relatively well, and no augmented sample is needed in the ALEC emulator. For the other three functional groups, the active learning algorithm has the smallest predictive RMSE on density ρ and grand potential Ω among all four designs. It is worth noting that the number of training inputs from active learning is always not larger than that in the RS2 design, but the predictive error by active learning is much smaller than the one by RS2. This is because the active learning strategy accurately identifies test samples hard to be predicted, and a numerical solver is called for computing these samples. Thus the computing budget is more efficiently spent with respect to the accuracy of predictions. Figure 6 presents more detailed prediction results of the particle densities and grand potentials for the fourth class of external potentials. The boxplot in Fig. 6(a) gives the out-of-sample RMSE for predicting the particle densities. Many densities are inaccurately predicted based on the samples from RS1, RS2, or D-opt, resulting in large errors in predicting the grand potential ( Fig.  6(b)). In comparison, the ALEC emulator identifies the samples that are difficult to be predicted and it uses them to enhance the emulator, thus providing accurate predictions for the remaining test samples. Although the total computational cost of the ALEC is similar to the methods with RS2 and D-opt schemes, the ALEC emulator results in better predictive accuracy.

B. Predicting particle densities and grand potentials for all groups of external potentials
In this subsection, we describe a more challenging scenario, where we aim to predict four groups of samples together. We assume that only the chemical potential and the external potential are used as input, while the parametric forms of the external potential are not known. The density profile and grand potential are harder to be predicted accurately in this case, as the input contains more degrees of freedom, leading to larger variability of the density and energy functionals, compared to the ones from any single class of input function alone. Since the inputs are a random sample from four classes of external potentials, we use the first n ini = 80 samples as initial inputs to build the ALEC emulator. When a test sample comes, we use the criterion from Eq. (15) to determine whether the test sample can be accurately predicted by the ALEC emulator or a numerical solver will be needed for solving the system directly.
We show the emulator's overall performance in the leftmost group (named as all groups) in Fig. 7 (detailed results are displayed in the supplementary materials). On average, the ALEC emulator performs much better, where the predictive RMSE of both particle densities and grand potentials is at least one magnitude smaller than all other approaches. A total of 873 training samples are added through the active learning algorithm, which is more than the total number of training samples added for building the emulator separately. This is because having a larger input space containing all 4 classes of external potentials increases the difficulty of learning densities, requiring more samples given the same threshold of the error. Note that the good predictive performance by ALEC is achieved using a similar or smaller sample size in training the emulator compared to the D-opt and RS2 approaches, by efficiently allocating the computing budget on computing particle densities hard to be predicted, FIG. 7. Out-of-sample RMSE of particle density ρ (bar plots) and grand potential Ω (line plots) for active learning with error control, RS1 design, RS2 design and active learning with D-optimality. The overall RMSE for combined all groups (solid bars and dots) and the RMSE of each individual group (shadow bars and dots) are both recorded. ALEC contains 80 inputs as initial training size, and the total number of augmented samples in the ALEC emulator (δ = 0.01, α = 0.05) is 873, with 0, 1, 58, and 814 augmented samples in groups 1 to 4, respectively. In D-opt, the total number of augmented samples is 935, and the augmented sizes for groups 1-4 are 59, 73, 393, and 410, respectively. RS1 contains a total of the first 80 training samples after uniformly mixing all four groups of density together. RS2 contains the first 953 training inputs, with 245, 221, 238, and 249 for groups 1-4, respectively. Since the ALEC emulator identifies the hardest region (group 4) to be predicted, thereby augmenting most samples from group 4, it has the smallest predictive error overall and predictive errors of all groups are controlled. and accurately predicting the rest of the densities with a negligible computational cost.
When examining the performance of different approaches in each function group, we found that the RMSE of particle densities from ALEC is homogenous across groups, while RS2 and D-opt approaches have smaller predictive errors for the first two groups, but have much worse performance for the fourth group. This is because RS2 and D-opt select more points from the first two groups to train the emulator, thereby yielding a similar or better result than ALEC in these two groups. However, the error from the first two groups can be controlled below the threshold of 0.01 for more than 95% of the test samples with no augmented samples or only a small number of augmented samples. In the ALEC emulator, for instance, the augmented samples from groups 1 and 2 are 0 and 1, respectively. On the contrary, the particle densities and grand potential energy from group 4 are the hardest to be predicted. This is automatically identified by ALEC, and the most computing budget (814 augmented samples) is spent on the fourth group to ensure the predictive errors of most samples are controlled below the threshold. This is exactly what we intend to accomplish. The ALEC can identify input regions that are difficult to predict, and then put more computing resources on those regions to control the predictive error. Consequently, the predictive error is controlled below the threshold for all subclasses, even if we do not know the labels of subclasses from the test samples, as all input classes are combined for predictions.
Another interesting finding is that the additional training samples from other input functional classes do not substantially improve the predictive performance for a specific input class. This is because the inputs from one class are distant from another class, as shown in Fig. 1(b), and thus the correlation between different input classes is small. It may be of interest to automatically identify the cluster in the training and split the training sample when constructing the emulator. This is useful as splitting the samples can reduce the computational complexity of the emulator. To explore this idea, we utilized a simple input decomposition technique. When the size of a cluster reaches 400, we divide it into two sub-clusters using the K-nearest neighbors (KNN) algorithm 39 . In this approach, the number of augmented samples and predictive RMSE of particle densities are 779 and 9.73 × 10 −4 , respectively. It achieves similar accuracy with 94 fewer training samples than the ALEC algorithm without input decomposition. Developing a more comprehensive way for input decomposition approach can further enhance the efficiency and accuracy of the surrogate model. Fig. 8 compares the computational time required in the ALEC approach and in calling a numerical solver (NS) each time. The largest computational cost in the ALEC algorithm comes from calling the NS for computing the particle densities of the initial and augmented training samples, while the time in constructing the ALEC emulator and making predictions is negligible, as the training sample size is not large. Since the ALEC emulator only requires a fraction of the samples to be numerically solved, it is much faster than running a numerical solver for each test sample. Note that RS2 and D-opt compared herein have a similar computational cost as the ALEC emulator, as the number of density profiles to be solved by NS is similar, yet the predictive error was not controlled in those approaches.

C. Predicting a new class of external potential
In this section, we examine the performance of the ALEC approach to predict inputs from a new functional form that is not in the training data set. To test this case, in addition to the existing 8000 samples, we generated 2000 new samples based on the weighted inputs between groups 2 and 3 as explained in Appendix B. We test the performance of the ALEC emulator for predicting the new input class, starting with the final GP emulator developed for groups 1-4 and all groups combined.
The left panel of Fig. 9 displays empirical CDF plots of the absolute errors of the ALEC emulator using the average variance selection criterion in Eq. (15). First, the predictive error of most of the test samples in the new group is small, for any initial emulator built upon other groups. Less than α = 5% of the absolute errors are greater than the chosen error bound δ = 0.01 for models using groups 1-3 and all groups combined. When we start from the initial emulator built from group 4 input, 11.2% of predictive errors are greater than δ = 0.01, which is larger than the probability tolerance threshold α = 5%. A close examination reveals that the percentage of the absolute error greater than δ is typically large for coordinates at the tail of the particle densities ( Fig.  9 middle panel), because of a rapid change in density at both ends, making it distinct from existing densities in group 4. The difference in the particle densities in the new class and group 4 increases the model discrepancy and makes the uncertainty in predictions harder to be quantified. Among the 11% of the test samples with RMSE larger than the threshold, however, most RMSEs are still very close to the threshold value, making this criterion ideal if the computational budget is not large.
On the other hand, the right panel of Fig. 9 shows the CDF of the RMSEs based on the maximum variance selection criterion in Eq. (12). Much less than α probability of the predictive error exceeds the threshold δ for all the cases shown here, as the maximum variance selection criterion is conservative, yielding more augmented samples and higher predictive accuracy. When the model discrepancy is large, one may use the maximum variance selection criterion to achieve a higher level of accuracy in predictions.

V. CONCLUSION
In this work, we propose the ALEC emulator, an active learning framework for controlling the predictive error of approximating computationally expensive functionals, based on the internal uncertainty assessment of the statistical surrogate model. This result has a probabilistic interpretation and it is fully automatic. We test the ALEC emulator on cDFT calculations of one-dimensional hard rods with different chemical potentials and external potentials. Through the numerical simulations of the density profiles and grand potentials for different groups of functions, we show that ALEC outperforms the conventional random sample designs, as well as the active learning approach with the D-optimality criterion in terms of prediction accuracy. With a fraction of the computational cost of running a numerical solver for each test input, we are able to predict functionals accurately with a controlled predictive error using the ALEC emulator. The ALEC emulator is designed as a reliable building block to be integrated for solving challenging simulation tasks, such as predicting expensive cDFT calculations on functional input, since the predictive error of the ALEC emulator is controlled for the test samples.
There are several future research directions. First, the ALEC emulator was built upon the PP-GP surrogate model for predicting vectorized output, while the approach can be generalized to using other surrogate models with uncertainty assessment in predictions. Physical symmetries, for instance, can be incorporated to define a new descriptor and distance metric to constrain the surrogate model. For instance, a common way to ensure translational invariance is to use pairwise distances of input, instead of the input function itself. Existing functionals found in physics 40 may be integrated as the mean of the emulator to improve its accuracy. Furthermore, dimension reduction of the input space through orthogonal representation can be helpful for reducing computational costs and improving the accuracy of predicting systems on 3D coordinates. Potential approaches include basis representation for reducing the dimensionality of the output space 30 and orthogonal projection of the inputs by maintaining large gradients of outcomes for reducing the dimensionality of the input space 41 . Finally, one may jointly model the grand potential energy and densities for predicting both quantities. This is particularly useful for applications with complex molecular systems, where computing the potential energies from particle densities is expensive.

SUPPLEMENTARY MATERIALS
The supplementary materials contain the derivation of predictive t-distribution in Eq. (5), and additional prediction results of particle density and grand potential. where a is the length of each hard rod.

Appendix B: External potentials and chemical potentials
Let L be the distance between two hard walls, and a is the length of the hard rods (HR). We consider the following four parametric forms of the external potential in this paper.
In the second functional form, the hard rods experience a van der Waals-like attraction from the hard walls (B2) Eq.(B2) is adopted from 43 . With fixed a and L, is a parameter in the range of ∈ (0.1, 2.2). Note that when = 0, we have the same potential as the first one. The third functional form consists of hard-wall potentials and a repulsion linearly proportional to the distance V ext (s) = ∞, s < a/2, s > L − a/2, m g s, a/2 < s < L − a/2.
Intuitively, the linear form mimics the gravitational potential. 44 This is the only asymmetric potential considered in this work. In training the machine-learning models, the range of m g is chosen to be (0.1, 3).
The power-law external potential 44 is the fourth functional class considereed in this work. This external potential has more flexibility as it has three parameters, u 0 , x 0 , and a 0 that can be used to modify the range of interactions: V ext (s) = ∞, s < a/2, s > L − a/2, u 0 s−L/2 x0 a0 , a/2 < s < L − a/2, where a 0 > 0. In this work, we use u 0 ∈ (1, 3), x 0 ∈ (1, 3) and a 0 ∈ (2, 5). Fig. 10 shows five examples of generated densities for each external potential.