Eigenvector Continuation as an Efficient and Accurate Emulator for Uncertainty Quantification

First principles calculations of atomic nuclei based on microscopic nuclear forces derived from chiral effective field theory (EFT) have blossomed in the past years. A key element of such ab initio studies is the understanding and quantification of systematic and statistical errors arising from the omission of higher-order terms in the chiral expansion as well as the model calibration. While there has been significant progress in analyzing theoretical uncertainties for nucleon-nucleon scattering observables, the generalization to multi-nucleon systems has not been feasible yet due to the high computational cost of evaluating observables for a large set of low-energy couplings. In this Letter we show that a new method called eigenvector continuation (EC) can be used for constructing an efficient and accurate emulator for nuclear many-body observables, thereby enabling uncertainty quantification in multi-nucleon systems. On the basis of ab initio calculations for the ground-state energy and radius in 4He, we demonstrate that EC is more accurate and efficient compared to established methods like Gaussian processes.

Introduction In recent years significant progress has been achieved in the theoretical and algorithmic development of sophisticated many-body methods that allow to study atomic nuclei up to mass number A 100 (see, e.g., Refs. [1][2][3][4][5][6] and references therein) based on nucleonnucleon (NN) and three-nucleon (3N) interactions derived from chiral EFT [7][8][9]. Given these many-body advances, the development of novel and more accurate nuclear interactions is a very active field of research. In addition to the theoretical work towards understanding how nuclei emerge from EFTs of the strong interaction, much effort is spent on the calibration of model parameters, e.g., low-energy constants (LECs) in EFT descriptions of nuclear interactions. In principle, calculations based on such interactions allow for a rigorous quantification of theoretical uncertainties stemming both from the parameter-estimation procedure as well as from truncating the EFT expansion at a given order. A rigorous uncertainty analysis is certainly possible and requires a careful determination of relevant covariances [10][11][12] and subsequent error propagation in all model predictions. Recently, Bayesian inference has been identified as a powerful and versatile tool for statistical analysis of EFTs, see for example Refs. [13][14][15][16][17][18][19][20].
Both parameter estimation and the calculation of posterior probability distributions for nuclear EFT or model predictions typically require extensive numerical sampling in a high-dimensional parameter space. Except for the simple two-nucleon sector, repeated calculation of nuclear many-body observables quickly becomes prohibitively expensive to allow for sample sizes sufficiently large to be meaningful. Yet, there are clear indications that many-body observables contain useful information for calibrating nuclear forces. For example, a fit of LECs to nuclear data including binding energies and radii of selected oxygen and carbon isotopes [21] showed that exploiting the information content of complex observables is phenomenologically important. In a similar spirit, input from α-α scattering data has been used to constrain two-nucleon forces [22,23]. In addition, it is clear that at least three-nucleon forces are necessary for an accurate theoretical description of nuclear systems based on EFT interactions. The LECs that enter for multi-nucleon forces need to be determined using calculations of light nuclei (typically A = 3, 4 are used), and already such calculations can incur a significant computational cost when a large number of them is needed. This significant computational cost highlights the importance of developing fast and accurate methods that make it possible to sample large parameter spaces using emulators, i.e., calculations that sacrifice the accuracy of an exact calculation for a significant gain in speed. The simplest such method, polynomial interpolation between a set of points within the parameter space, is usually not a viable option for a lack of both accuracy and efficiency. Gaussian processes [24] are useful for leveraging expensive statistical analyses in nuclear theory [20,25]. As a machine-learning method they can be advantageous for systematically exploring large parameter spaces and by design provide uncertainties of the emulator output, but like polynomials they are still limited to interpolation within a set of training data and cannot be used for reliable extrapolations. In this Letter, we explore eigenvector continuation (EC), introduced in Ref. [26], as an alternative to overcome this limitation. We find that EC performs accurate extrapolations in multi-dimensional parameter domains even to points far outside the training data set used to construct the emulator, and that it pro- vides a significantly more efficient and accurate emulator of nuclear systems than a Gaussian process.
Formalism Eigenvector continuation is based on the fact that when a Hamiltonian depends smoothly on some real-valued control parameter, any eigenvector of the Hamiltonian is a smooth function of that parameter as well. Furthermore, the eigenvector trajectory traced out as the parameter is varied can be well approximated by a finite-dimensional manifold [26]. This last statement can be turned into a variational method for computing the eigenvector for any value of the control parameter.
Consider a Hamiltonian H(c) that varies smoothly with real parameter c.
The ground-state eigenvector |v 0 (c) can be well approximated as some linear combination of the ground-state eigenvectors |v 0 (c [1] ) , · · · , |v 0 (c [N ] ) at "training points" c [1] , · · · , c [N ] . In order to determine the desired linear combination that best approximates |v 0 (c) , we simply find the ground state of H(c) projected onto the subspace spanned by |v 0 (c [1] ) , · · · , |v 0 (c [N ] ) . In Ref. [26] the applications of EC focused mainly on extrapolation in cases where the direct calculation of |v 0 (c) was not possible due to computational issues such as the Monte Carlo sign problem. In this work we will use EC for both interpolation and extrapolation. We also consider, for the first time, the extension of EC to Hamiltonians that depend on more than one control parameter.
Specifically, we explicitly demonstrate the advantages of using EC for constructing a fast and accurate emulator for nonrelativistic calculations of the 4 He nucleus. While this application is a benchmark case that is particularly relevant for nuclear physics, the very general mathematical underpinnings of EC enable the emulation of expensive problems across several disciplines also outside of physics provided only that they can be formulated as an eigenvalue problem. Eigenvector continuation moreover supports full reconstruction of the emulated eigenvector (wavefunction). To demonstrate this we consider both the ground-state energy E and radius r of 4 He nucleus, as functions of the 16 LECs c in a particular chiral potential V (c) for the strong interaction [12], entering the Schrödinger equation H(c)|ψ(c) = E|ψ(c) .
Training the EC emulator consists of building a basis to span an eigenvector subspace. For this we must obtain exact eigenvectors (wavefunctions) |ψ(c [i] ) for a set of N EC points c 1 , . . . , c NEC across the chosen 16dimensional parameter domain of the LECs. We formulate the Schrödinger equation for 4 He as an eigenvalue equation using the no-core shell model (NCSM) [27]. This is a variational basis-expansion method, also known as "configuration interaction" in quantum chemistry. The exact wave function of the Hamiltonian H(c i ) is expanded in eigenfunctions of the harmonic-oscillator (HO) potential, yielding a Hamiltonian represented as a matrix in this HO basis that is subsequently diagonalized. Considering low-energy states motivates a truncation of this expansion based on a maximum number of oscillator quanta N max . Another parameter characterizing the basis is the oscillator frequency Ω. For N max → ∞, the choice of frequency is arbitrary, but for each truncated basis there is a residual dependence of results on Ω that has to be assessed [28,29]. The underlying many-body problem is translationally invariant and thus preferably expressed in relative coordinates. For few-body systems like 4 He it is possible to proceed this way, which includes an exact evaluation of the four-fermion antisymmetrizer. For systems with more than four nucleons, it is however computationally more efficient to antisymmetrize in single-particle coordinates [30]. To leverage a comparison between the EC emulator and exact solutions we truncate the HO basis expansion at N max = 16 for a frequency Ω = 36 MeV, which typically gives sub-percent accuracy for the ground-state energy and radius of 4 He. With this choice the HO basis consists of 2775 antisymmetric and translationally invariant four-body states.
The nuclear potential that we employ is additive in the d = 16 LECs, i.e., we can express the Hamiltonian as where H 0 includes the kinetic energy. Any Hamiltonian with more than one interaction parameter can be written in this form. Furthermore, each term H i for i = 1, . . . , 16 can be projected onto the EC subspace once and then used for an arbitrary number of emulations. Each of these corresponds to a straightforwards solution of the N EC ×N EC -dimensional generalized eigenvalue problem.
Results To systematically investigate the quality of the EC emulator, we consider several different cases for the number of LECs that we vary simultaneously, amounting to sampling Hamiltonians in a d-dimensional parameter space, where d = 1, . . . , 16. We select the set of training points T = {c [i] } NEC i=1 using a space-filling Latin Hypercube design [31]. For simplicity we define a parameter domain for each LEC between −2 and 2 in appropriate units of inverse energy, see, e.g., Ref. [21]. Validation data is drawn randomly from a uniform distribution U(−2, +2). Each validation point c corresponds to either interpolation or extrapolation from the set of training points, with the former being defined as the case where c lies within the convex hull of T . By randomly generating a coefficient vector α with α k ≥ 0 for k = 1, . . . , d and k α k = 1 it is possible to alternatively sample only points k α k c [k] corresponding to interpolation. We present results as a cross-validation plots where   Fig. 2, i.e., varying 16 LECs and using an EC subspace constructed from 64 training data points. The assumed number of matrix-vector products required for a Lanczos diagonalization in the full Nmax = 16 space is Nmv = 80 for this case (see Appendix and main text for details). The theoretical limit indicates the max speedup reached asymptotically in the number of samples, which is 614 in the present case.
we consider emulated values as a function of the exact ones. In these plots we include results for polynomial interpolation and a Gaussian process for comparison. The Gaussian process is constructed using a standard squared exponential kernel with hyperparameters estimated from the maximum of the marginal log-likelihood of the calibration data. A Python script able to run calculations of this type is provided along with this Letter, with a brief explanation included in the Supplemental Material.
A representative example is shown in Fig. 1. In this case, calculations for the 4 He ground-state energy are emulated as a function of three LECs using 12 training data points obtained in an N max = 16, Ω = 36 MeV NSCM model space. Eigenvector continuation is seen to work exceptionally well (the difference to exact calculations for each point is negligibly small and cannot be resolved in the plot), whereas polynomial interpolation and the Gaussian process struggle to provide accurate results even when we consider only validation points corresponding to interpolation within the convex hull of the set of training points (right panel in Fig. 1).
In fact, EC can achieve excellent results even with fewer than 12 training data points in this particular case. Furthermore, EC requires only a moderate increase in the number of training data as the dimension of the parameter space is increased. In Fig. 2 we show results for the 4 He energy with all 16 LECs varied, using the same N max = 16, Ω = 36 MeV NSCM model space as before. It is evident how EC can still provide accurate results while polynomial interpolation and the Gaussian process  fail completely to emulate the data, even though only interpolation is considered in Fig. 2.
To fully appreciate the efficiency gain provided by the EC method, it is important to compare the overall computational cost of the different methods considered above. The cost of emulating with EC is not severe because all relevant matrix operations, i.e., setting up the target Hamiltonian and solving a generalized eigenvalue problem, need only be performed in the small EC subspace. Besides the requirement of carrying out N EC exact calculations there is a one-time cost of matrix-matrix-matrix multiplications coming from projecting the Hamiltonian to the EC subspace. Thus, the benefit of emulating with EC will improve with the number of calls to the emulator. Asymptotically in the number of emulator calls, the speedup of using EC is proportional to (M/N EC ) 2 , where M is the dimensionality of the full-space problem. Typically, we find N EC ≈ 10 − 100 for problems with M ≈ 10000, thus easily yielding a speedup factor ∼ 10 4 . In Fig. 3 we show the speedup we achieved for the 4 He problem benchmarked here. A detailed analysis of the computational cost is provided in the Supplemental Material.
Gaussian-process interpolation provides an uncertainty estimate of the output value by design. For EC, a detailed analysis of its rate of convergence as the number of training points is increased is work in progress [32]. However, even without a fully developed theory for EC uncertainties, we can make the following remarks: First, EC is a variational method. This can be seen directly by noting that it is based on constructing a subspace: considering the original Hamiltonian in diagonalized form, it is clear that removing any of the basis vectors can only increase the lowest eigenvalue of the remaining operator. Therefore, EC-emulated (energy) eigenvalues will always be larger than or equal to the true result, i.e., resulting in one-sided error bars. Note, however, that this argument does not apply to other operators evaluated in the EC subspace. Second, the fact that EC provides remarkably accurate results with only a small amount of training data, as well as the benefit that it can reliably both interpolate and extrapolate, can be exploited in order to estimate the uncertainty based on removing different points from the training set, giving a range of values for each emulation target point. For many applications, this may be an efficient strategy to assess how converged the EC-emulated values are. For a more thorough analysis one can use various training sets of the same size and analyze the distribution of results. We have used this strategy to obtain the results for the 4 He ground-state radius squared shown in Fig. 4. Specifically, the figure shows the uncertainty bands obtained by considering 32 additional training data sets of 128 points each. The thicker uncertainty bars correspond to the results obtained from 68.2 % of all sets, while the faint thinner ones indicate the full range of results for each point. The results indicate that the size of the resulting uncertainty bars correlate well with the degree of deviation from the exact results and hence serve as a possible reasonable estimate for the uncertainties.
Conclusion and outlook We have demonstrated how EC can be used to construct an efficient and accurate emulator of eigenvalue problems with continuous and highdimensional parametric dependencies. Moreover, for systems with a matrix representation that linearly depends on a set of parameters, the EC method enables a substantial computational speedup while maintaining highaccuracy outputs compared to exact solutions of the original problem. This is achieved by constructing a tailored low-dimensional subspace spanned by exact eigenvectors for a set of "training" points in the parameter space. We constructed an efficient and accurate emulator of the quantum-mechanical solution of the 4 He nucleus, considering its ground-state and squared radius as concrete observables. The computational speedup offered by the EC emulator is essential for sampling high-dimensional regions in the parameter domain of any model with the purpose of, e.g., optimization and uncertainty quantification, where the required large number of exact calculations would be prohibitively expensive. For nuclear physics, the EC method can be a key ingredient to facilitate large-scale Markov-Chain Monte Carlo evaluations of relevant Bayesian posteriors of the parameters in EFTs or models of the nuclear forces. Applications to this and related studies are already under way.

Cost comparison
For the following analysis, we let M = M (N max ) denote the actual dimension of the model space considered in a given calculation (suppressing the dependence on the number of nucleons). Furthermore, N EC is the number of training data points, i.e., the number of states spanning the EC subspace, while N denotes the number of requested samples. We assume that sufficient memory is available to store intermediate results as necessary and we limit the analysis to basic estimates for the required operations, not taking into account specific optimizations that may be used in practice.
• We first consider the cost of performing a single calculation in the full M -dimensional space. Setting up the Hamiltonian, given by a part independent of LEcs plus a linear combination of terms for each individual LEC, H = H 0 + NLEC α=1 c α H α , costs a total of 2N LEC M 2 floating-point operations. Subsequently calculating the ground-state energy with a Lanczos-like algorithm has a complexity that is dominated by performing N mv Mdimensional matrix-vector multiplications, each of which costs M 2 operations. Note that the specific value of N mv depends on the desired accuracy of the calculation as well as on the properties of the Hamiltonian. In particular, N mv typically grows with increasing M . Neglecting other aspects of the diagonalization procedure, we arrive at a total cost of M 2 × (2N LEC + N mv ) operations.
• Multiplying the above by N gives the cost for a direct sampling within the full space.
• Setting up an emulator has a base cost of N EC × M 2 × (2N LEC + N mv ). For polynomial interpolation and Gaussian process emulation we take this as the total cost and assume the subsequent cost for obtaining samples is negligible.
• Setting up sampling based on EC requires some additional work. • For each point sampled using EC, the Hamiltonian setup then only needs to be performed in the subspace, amounting to 2N LEC N 2 EC operations per sample. Solving the generalized eigenvalue problem costs another 14N 3 EC operations [33].
• The sampling cost can be reduced by performing an initial orthogonalization of the {c i } NEC i=1 (which we assume to be achieved through a singular-value decomposition costing about 6M N 2 EC + 11N 3 EC operations [33]), leaving only the solution of a standard symmetric eigenvalue problem and thus a cost of 26N 3 EC /3 + O(N 2 EC ) operations per sample [34].

Python code
We provide the Python code ec xval.py as supplementary material. This program is a simplified version of the script that was used to generate the cross-validation plots shown in the main text. Matrices required as input data (NCSM Hamiltonian along with corresponding representation of radius squared operator) are provided as well. Due to storage limitations, these matrices are restricted to rather small NCSM model spaces, but they nevertheless provide representative examples. It is our hope that this code will facilitate applications of eigenvector continuation to a variety of cases where efficient and accurate emulators are required. Making use of freely available Python packages, the code generates cross validation plots that compare EC to both a Gaussian process and simple polynomial interpolation.
Running a cross validation is as simple as typing $ python3 ex_xval.py -d 3 -n 6 in the terminal. This will generate a cross-validation plot for a three-dimensional parameter space, using 6 EC basis vectors. By default, the cross validation is run using only eigenvector continuation. In order to compare at the same time to polynomial interpolation and a Gaussian process, as shown in the main text, -pge can be given as an option to enable all emulators. A number of further aspects can be controlled by passing command-line options, a full list of which, along with explanations, is printed to the terminal by running: