# Training and Projecting: A Reduced Basis
Method Emulator for Many-Body Physics

Edgard Bonilla  Bonilla and Giuliani are co-authors contributing equally to
this work<EMAIL_ADDRESS>Department of Physics, Stanford University,
Stanford, California 94305, USA  Pablo Giuliani<EMAIL_ADDRESS>FRIB/NSCL Laboratory, Michigan State University, East Lansing, Michigan 48824,
USA  Department of Statistics and Probability, Michigan State University, East
Lansing, Michigan 48824, USA  Kyle Godbey<EMAIL_ADDRESS>FRIB/NSCL
Laboratory, Michigan State University, East Lansing, Michigan 48824, USA  Dean
Lee<EMAIL_ADDRESS>FRIB/NSCL Laboratory, Michigan State University, East
Lansing, Michigan 48824, USA  Department of Physics and Astronomy, Michigan
State University, East Lansing, Michigan 48824, USA

###### Abstract

We present the reduced basis method (RBM) as a tool for developing emulators
for equations with tunable parameters within the context of the nuclear many-
body problem. The RBM uses a basis expansion informed by a set of solutions
for a few values of the model parameters and then projects the equations over
a well-chosen low-dimensional subspace. We connect some of the results in the
eigenvector continuation literature to the formalism of RBMs and show how RBMs
can be applied to a broad set of problems. As we illustrate, the possible
success of the RBM on such problems can be diagnosed beforehand by a principal
component analysis. We apply the RBM to the one-dimensional Gross-Pitaevskii
equation with a harmonic trapping potential and to nuclear density functional
theory for 48Ca. The outstanding performance of the approach, together with
its straightforward implementation, show promise for its application to the
emulation of computationally demanding calculations, including uncertainty
quantification.

Most modern theoretical models describing many-body nuclear dynamics share an
ever-increasing computational burden. This can turn into a challenge tasks
like uncertainty quantification (UQ) analysis Phillips  _et al._ (2021);
Ekström  _et al._ (2019), experimental design Melendez  _et al._ (2021a);
Giuliani and Piekarewicz (2021), calibration of model parameters Pratt  _et
al._ (2015); McDonnell  _et al._ (2015); Wesolowski  _et al._ (2019), and even
repeated evaluation for different inputs Lasseri  _et al._ (2020).
Emulators—algorithms able to give fast-and-accurate approximate calculations
to expensive computations—have been gaining increasing importance as a way to
circumvent these challenges Santner  _et al._ (2003); Phillips  _et al._
(2021).

In recent years, a technique called eigenvector continuation (EC) Frame  _et
al._ (2018) was developed to emulate computationally intensive calculations
involving bound states of Hamiltonian operators König  _et al._ (2020) and
nuclear scattering Furnstahl  _et al._ (2020); Melendez  _et al._ (2021b);
Drischler  _et al._ (2021). EC has shown excellent performance in
interpolation and extrapolation by working with two elements: choosing its
basis functions (the ansatz) from the linear span of exact solutions to the
problem at hand, and using a variational principle—for example the Rayleigh-
Ritz method Frame  _et al._ (2018); Cohen-Tannoudji  _et al._ (1986a) or the
Kohn variational principle Furnstahl  _et al._ (2020); Kohn (1948)—to obtain
equations for the coefficients of this linear combination.

We present an emulator constructed in the formalism of reduced basis methods
(RBMs) Almroth  _et al._ (1978); Quarteroni  _et al._ (2015); Hesthaven  _et
al._ (2016), a set of dimensionality-reduction techniques that have proven
useful in a variety of computationally-intensive problems involving partial
differential equations Quarteroni  _et al._ (2011); Nguyen  _et al._ (2010);
Field  _et al._ (2011); Milani  _et al._ (2008). EC can be naturally connected
to RBMs by constructing a generalization of EC through a Galerkin method
formulation. The key insight is that once a reasonable choice of ansatz
functions has been made (for example, the EC basis), all that is needed is a
method to select a suitable candidate approximation from the ansatz subspace.
This could be achieved, for instance, by using a variational principle,
minimizing a cost functional, or by finding the fixed point of an iterative
scheme. Among the alternatives, the Galerkin method—the option often chosen in
RBMs—stands out for its simplicity: it attempts to find an accurate
approximate solution by projecting the problem to a well-chosen very-low
dimensional subspace. This simplicity allows the RBM formulation to be applied
to a wide variety of problems in a straightforward way.

RBMs are tailored to problems that feature an equation that depends smoothly
on a list of tunable control parameters $\alpha$Quarteroni  _et al._ (2015).
The goal is to build an approximate solution for a suitable range of these
parameters. Let us assume the equation is written in the general form:

$F_{\alpha}(\phi_{\alpha})=0,$ (1)

where $\phi_{\alpha}$ is a vector (or function) from a Hilbert space
$\mathcal{H}$, and $F_{\alpha}$ maps $\mathcal{H}$ onto itself. For example,
in the case of bound systems with a Hamiltonian $H_{\alpha}$ that depends on
$\alpha$, $F_{\alpha}$ can take the form of the eigenvalue equation
$F_{\alpha}(\phi_{\alpha})=H_{\alpha}\phi_{\alpha}-\lambda_{\alpha}\phi_{\alpha}$,
where $\lambda_{\alpha}$ is the eigenvalue. Another example would be the case
of single-channel scattering where $F_{\alpha}$ can be the radial part of the
scattering equation Thompson and Nunes (2009)
$F_{\alpha}(\phi_{\alpha})=\left(-\frac{d^{2}}{dr^{2}}+\frac{\ell(\ell+1)}{r^{2}}+U(r,\alpha)-p^{2}\right)\phi_{\alpha}(r)$,
where a system with reduced mass $\mu$ interacts through a potential
$V(r,\alpha)=U(r,\alpha)/2\mu$ with parameters $\alpha$, $\ell$ is the angular
momentum quantum number, and $p$ is the asymptotic linear momentum. The RBM
finds approximate solutions $\hat{\phi}_{\alpha}$ to these—and more
general—problems by constructing a basis expansion with $n$ linearly
independent ‘reduced basis’ (RB) functions $\\{\phi_{k}\\}_{k=1}^{n}$:

$\hat{\phi}_{\alpha}=\phi_{0}+\sum_{k=1}^{n}a_{k}\phi_{k},$ (2)

where $\phi_{0}$ is an extra term that can be added to satisfy boundary
conditions imposed on Eq. (1). The RB functions $\phi_{k}$ are selected to
create an affine space (the ansatz subspace) _close_ to the manifold formed by
the solutions to Eq. (1) as a function of the parameters $\alpha$ Quarteroni
_et al._ (2015) by using the information from a (possibly small) sample of
exact solutions. In practice, the ‘exact solutions’ (or ‘snapshots’) of Eq.
(1) are constructed by highly accurate yet computationally expensive
approximations such as finite element or spectral calculations Quarteroni  _et
al._ (2011).

An option for building the RB in Eq. (2) is the approach taken in many EC
applications Frame  _et al._ (2018); König  _et al._ (2020); Furnstahl  _et
al._ (2020); Melendez  _et al._ (2021b); Drischler  _et al._ (2021), also
known as the Lagrange RB Quarteroni  _et al._ (2011). It consists of
calculating $n$ ‘training functions’ $\\{\tilde{\phi}_{k}\\}^{n}_{k=1}$ as
solutions to Eq. (1) for $n$ values $\alpha_{k}$
($F_{\alpha_{k}}(\tilde{\phi}{{}_{k}})=0$), and then choosing the RB as these
$n$ training functions $\phi_{k}=\tilde{\phi}_{k}$. One possible way to
improve upon this choice is the proper orthogonal decomposition (POD)
Quarteroni  _et al._ (2015). It consists of computing $N\geq n$ solutions
$\\{\tilde{\phi}_{l}\\}^{N}_{l=1}$ and constructing the RB with the first $n$
components from a principal component analysis (PCA) Jolliffe (2002), or
singular value decomposition (SVD) Blum  _et al._ (2020), of the set of these
$N$ training functions. Therefore, by using the information of $N$ samples,
the POD basis is more robust than a Lagrange RB of dimension $n$, and faster
than a Lagrange RB of $N$ training points. As a side note, an exponential
decay on the associated singular values $\sigma_{k}$ from the PCA indicates
that the RBM can provide an accurate approximation for the problem at hand
Quarteroni  _et al._ (2011, 2015); Nonino  _et al._ (2019). We exploit this
feature later when discussing Fig. 1.

Once the $n$ RB functions are chosen, the coefficients $a_{k}$ for the
approximation are found by the Galerkin method Rawitscher  _et al._ (2018),
that is, by projecting Eq. (1) over $n$ linearly independent ‘projecting
functions’ $\\{\psi_{k}\\}_{k=1}^{n}$ in the Hilbert space:

$\langle\psi_{k}|F_{\alpha}(\hat{\phi}_{\alpha})\rangle=0,\quad\text{for all
}k.$ (3)

$F_{\alpha}(\hat{\phi}_{\alpha})$ is often called the residual Fletcher
(1984), and it can be used, for example, to inform the construction of the RB
Buffa  _et al._ (2012), or to estimate the emulation error Prud’homme  _et
al._ (2001). We can interpret Eq. (3) as enforcing the orthogonality of
$F_{\alpha}(\hat{\phi}_{\alpha})$ to the subspace spanned by
$\\{\psi_{k}\\}_{k=1}^{n}$, i.e., by finding a $\hat{\phi}_{\alpha}$ such that
$F_{\alpha}(\hat{\phi}_{\alpha})$ is “zero” up to the ability of the set
$\\{\psi_{k}\\}_{k=1}^{n}$. The choice of projecting functions $\psi_{k}$ is
arbitrary, but is usually also informed by the solution manifold Quarteroni
_et al._ (2015); Hesthaven  _et al._ (2016). For the rest of this work, we
choose $\psi_{k}$ to enforce orthogonality with respect to the ansatz subspace
(2), which is the traditional way of using the Galerkin method Fletcher
(1984).

To illustrate the application of this method and connect with previous results
in the EC literature, we work with the two previously mentioned examples for
$F_{\alpha}$. For the single-channel scattering example, we create the
approximate solution $\hat{\phi}_{\alpha}=\sum_{k=1}^{n}a_{k}\phi_{k}$, with
$\phi_{k}$ as exact solutions to $n$ different $\alpha$, and satisfy the
boundary conditions by imposing $\sum_{k=1}^{n}a_{k}=1$, as done in Furnstahl
_et al._ (2020). By letting $a_{1}=1-\sum_{k=2}^{n}a_{k}$ and thus
$\hat{\phi}_{\alpha}=\phi_{1}+\sum_{k=2}^{n}a_{k}(\phi_{k}-\phi_{1})$, we
explicitly satisfy the boundary condition while having only $n-1$ free
coefficients. We identify $(\phi_{k}-\phi_{1})$ for $k\geq 2$ as the relevant
elements in the basis expansion and select $\psi_{k}=(\phi_{k}-\phi_{1})$ as
the associated projecting functions, leading to the equations:

$\sum_{j=1}^{n}a_{k}\langle\phi_{k}-\phi_{1}|F_{\alpha}({\phi}_{j})\rangle=0,\quad\text{for}\quad
2\leq k\leq n.$ (4)

These equations—formulated from a geometric projection argument—are equivalent
to those obtained by the generalized Kohn variational principle Lucchese
(1989), as done in Drischler  _et al._ (2021). The proof is elaborated in the
supplementary material.

In the bound system example, with $\phi_{k}=\psi_{k}$, Eq. (3) might not have
a solution for the exact eigenvalue. Allowing $\lambda_{\alpha}$ to be
approximated by $\hat{\lambda}_{\alpha}$ helps ensure we can solve the
projected equations:

$\sum_{j=1}^{n}a_{j}\langle\phi_{k}|H_{\alpha}|{\phi_{j}}\rangle=\hat{\lambda}_{\alpha}\sum_{j=1}^{n}a_{j}\langle\phi_{k}|{\phi}_{j}\rangle,\quad\text{for
all }k,$ (5)

where the set $a_{k}$ plus the approximate eigenvalue $\hat{\lambda}_{\alpha}$
add up to $(n+1)$ unknowns. We complete the set of equations with a
normalization condition:
$\langle\hat{\phi}_{\alpha}|\hat{\phi}_{\alpha}\rangle=1$. When choosing
$\phi_{k}$ as exact solutions for different $\alpha$, Eq. (5) is equivalent to
the generalized eigensystem of EC formulated from the Rayleigh-Ritz
variational method Frame  _et al._ (2018).

Beyond these two examples, the generality of the Galerkin formalism allows to
apply RBMs to a wide variety of problems, including discrete, operator,
integral, and differential equations Fletcher (1984); Singh and Stauffer
(1974). As such, the projected equations (3) can be directly applied to non-
linear problems like non-linear eigenvalue equations where
$F_{\alpha}=G_{\alpha}({\phi}_{\alpha})-\lambda_{\alpha}\phi_{\alpha}$, with
$G_{\alpha}$ a general operator. Additionally, in the case of coupled
equations—common in many-body physics—the formalism is easily extended by
expanding each coupled function $\hat{\phi}^{(i)}_{\alpha}$ independently and
finding the coefficients $a_{k}^{(i)}$ by enforcing Eq. (3) with projecting
functions $\\{\psi_{k}^{(i)}\\}_{k=1}^{n}$ for each $i$-th coupled function.

It is important to note, however, that if the solution manifold
$\phi_{\alpha}$ for the problem at hand cannot be sufficiently embedded in a
linear subspace, then the RBM we described will not constitute an effective
emulator. We can engineer a simple example with the 1D quantum Hamiltonian of
a particle trapped in an infinite well (IW) Cohen-Tannoudji  _et al._ (1986b)
by letting $\alpha$ control the location of the well. A direct application of
the RBM fails to accurately emulate the ground state wave function as $\alpha$
changes. Extensions to the basic RBM can tackle these issues by allowing
further manipulation of the RB functions Nonino  _et al._ (2019).

In practice, to test whether a problem is fit for emulation via RBMs, it is
sufficient to observe the decay of the singular values $\sigma_{k}$ associated
with the PCA of a group of exact solutions $\tilde{\phi}_{\alpha}$ for various
$\alpha$ Quarteroni  _et al._ (2015). Fig. 1 shows the singular values
$\sigma_{k}$ for the problems discussed in our work, all decaying
exponentially except for the IW. In the case of the scattering wave functions
in the ${}^{1}S_{0}$ channel at a fixed energy for the Minnesota potential
Thompson  _et al._ (1977), the $\sigma_{k}$ are consistent with the EC results
Furnstahl  _et al._ (2020). A similar pattern is obtained for wave functions
across energies (black squares in Fig. 1) by re-scaling the scattering
differential equation via the change of variable $s=pr$, making all exact
solutions $\tilde{\phi}_{\alpha}(s)$ share the same asymptotic behavior. This
observation implies that in principle it should be possible to build a
scattering emulator across energies. For the Gross-Pitaevskii (GP) equation
and the 13 energy levels of 48Ca under density functional theory (DFT), the
decay of their respective $\sigma_{k}$ also makes them excellent candidates
for the application of the RBM, as we explore next.

Figure 1: Decay of the singular values $\sigma_{k}$ for a set of solutions of:
the Infinite Well (I.W. in magenta) with $V(x,\alpha)=0$ for
$x\in[\alpha,\alpha+1]$ and $V(x,\alpha)=\infty$ otherwise, the single channel
2-body Scattering with a Minnesota potential Thompson  _et al._ (1977) with
fixed energy (S.M. $E_{0}=50$ MeV in red) and varying energy (S.M.
$E^{*}\in[20,80]$ MeV in black) with the parameter’s range as used in Ref.
Furnstahl  _et al._ (2020), the Gross-Pitaevskii equation (G.P. in blue), and
the solutions for the 13 energy levels in 48Ca (DFT in dashed green-blue
lines). The supplementary material contains additional details on these
calculations.

The GP equation Gross (1961); Pitaevskii (1961) (see also Pichi  _et al._
(2020) for a RBM application) is a nonlinear Schrödinger equation that
approximately describes the low-energy properties of dilute Bose-Einstein
condensates. Using a self-consistent mean field approximation, the many-body
wave function is reduced to a description in terms of a single complex-valued
wave function $\phi(\vec{r})$. We work with the one-dimensional GP equation
Carr  _et al._ (2000a, b, 2001); Torres-Vega (2017) with a harmonic trapping
potential by letting $F_{\alpha}$ be:

$F_{q,\kappa}(\phi)=-\phi^{\prime\prime}+\kappa
x^{2}\phi+q|\phi|^{2}\phi-\lambda_{q,\kappa}\phi=0,$ (6)

where $\kappa$, $q$, and $\lambda_{q,\kappa}$ are proportional to the strength
of the harmonic trapping, the self-coupling of the wave function, and the
ground state energy, respectively. $\phi(x)$ is a single variable function
that depends on $x$ and it is normalized to unity.

Figure 2: Comparison between the exact solvers (solid lines) and the RBM
(points) calculations. Panel a) shows the ground state energy
$\lambda_{q,\kappa}$ of the GP equation as a function of $q$ for the values of
$\kappa=[0.5,2,4,10,30]$. Panel b) shows the neutron skin thickness of 48Ca as
a function of $L$ for the values of $K=[160,210,260]$ MeV. In both figures,
the Lagrange RB with four points is shown as cyan stars within the inset
plots. The supplementary material contains additional details on these
calculations.

To test the RBM for extrapolation, we built a Lagrange RB with four training
functions $\tilde{\phi}_{i}$ in the $[q,\kappa]$ space as exact solutions
($F_{q_{i},\kappa_{i}}(\tilde{\phi}_{i})=0$), with the projecting functions as
$\psi_{i}=\phi_{i}$. Panel a) in Fig. 2 shows the results of emulating
$\lambda_{q,\kappa}$ by using this RB and applying Eqs. (3) plus the
normalization condition. The agreement between the exact and RBM calculations
is excellent, with an error of less than $2.5\%$ in the repulsive phase
($q\geq 0$) where the four training parameters are located, and it
deteriorates only in the attractive phase ($q<0$) well beyond the training
region. Extrapolation is not a feature usually exploited on the RBM
literature, yet it could be key when calculations of exact solutions in a
specific phase of the system are numerically unstable or impossible, but
approximable by the RBM.

In addition to extrapolating, we explored a situation similar to how emulators
are tested for UQ König  _et al._ (2020); Drischler  _et al._ (2021). Using a
Latin hypercube sampling (LHS) Mckay  _et al._ (2000), we drew 500 testing
points in the range $q\in[0,30]$ and $\kappa\in[5,30]$. We constructed three
types of RB: Lagrange, POD, and POD+Greedy, each with three sizes $n=(2,4,8)$.
The Lagrange RB consisted of $n$ exact solutions drawn with LHS. The POD and
POD+Greedy consisted of $n$ principal components from a set of $N=20$ exact
solutions. For the POD the $N$ exact solutions were drawn using LHS, while for
the POD+Greedy the first solution was placed at a central location and the
other $N-1$ were included one-by-one through a Greedy algorithm inspired on
Refs. Veroy  _et al._ (2003); Haasdonk and Ohlberger (2008); Sarkar and Lee
(2021). Our Greedy approach finds the parameter set $[q_{m+1},\kappa_{m+1}]$
for the next exact solution $\tilde{\phi}_{m+1}$, by maximizing the norm of
the residual $F_{q,\kappa}(\hat{\phi}_{q,\kappa})$ over a LHS of parameters
$[q,\kappa]$. In each step, $\hat{\phi}_{q,\kappa}$ is constructed with a POD
basis informed by the previous $m$ exact solutions.

Table 1 shows the relative root mean squared errors (RMSEs), which converge
exponentially as expected from the results in Fig. 1. Both POD bases were more
accurate and robust than the Lagrange RB, which produced results that
frequently changed by more than an order of magnitude when re-sampling the
exact solutions for the basis. For $n=8$ the POD+Greedy algorithm’s accuracy
was more than $600$ times better than the Lagrange RB. In terms of speed-up
when calculating the 500 testing points, the three RB with $n=2$ were almost
$150$ times faster than the exact solver, while $n=4$ and $n=8$ obtained
speed-ups of 40 and $5$ times, respectively.

Table 1: RMSEs for the GP and DFT problems described in the text. The RMSE is defined as $\langle\big{[}(A_{\text{RBM}}-A_{\text{exact}})/A_{\text{exact}}\big{]}^{2}\rangle^{1/2}$, where $A$ is the quantity being computed, and $\langle\rangle$ denotes average. Three cases of the RB size were explored with $n=(2,4,8)$. 500 testing points were drawn in their respective parameter space, but for DFT 32 points were excluded from the statistics since the exact solver reported convergence problems. The supplementary material contains additional details on these calculations. Basis n |  | G.P. Ground
---
State Energy
|  48Ca Average
---
Particle Energy
Lagrange | POD |  | POD
---
Greedy
POD
2 | $1.0\times 10^{-1}$ | $1.2\times 10^{-2}$ | $1.5\times 10^{-2}$ | $5.9\times 10^{-3}$
4 | $3.0\times 10^{-3}$ | $5.6\times 10^{-4}$ | $2.1\times 10^{-4}$ | $6.1\times 10^{-4}$
8 | $1.3\times 10^{-5}$ | $1.2\times 10^{-6}$ | $2.0\times 10^{-8}$ | $1.7\times 10^{-4}$

We now proceed to use the RBM in realistic nuclear DFT calculations. DFT is a
widely applied microscopic formalism Schunck (2019) (see also Cancès  _et al._
(2007); Lin  _et al._ (2012); Zhang  _et al._ (2017) for other RBMs
applications to DFT). In nuclear physics it is used to describe properties of
nuclei from the mean-field perspective, i.e., each nucleon interacts with an
average effective field made up of all the particles in the system. This
interaction is then constructed in a self-consistent way: the wave function of
each nucleon and its eigenenergy are found at the same time as the effective
field they produce and interact with. As such, the Hamiltonian $\hat{h}^{(j)}$
acting on the $j$ wave function $\phi^{(j)}$ depends on all $M$ of them:

$\hat{h}^{(j)}[\Phi]\phi^{(j)}-\lambda^{(j)}\phi^{(j)}=0\quad\text{for}\quad
1\leq j\leq M,$ (7)

where $\Phi=\\{\phi^{(i)}\\}_{i=1}^{M}$, and the parameter list $\alpha$ has
been omitted for the sake of clarity. The dependence of the Hamiltonian on the
wave functions comes from, for example, the total nuclear density $\rho$ and
kinetic energy density $\tau$. We derive the single particle Hamiltonian,
$\hat{h}^{(j)}$, from the Skyrme effective interaction Skyrme (1956, 1958);
Vautherin and Brink (1972); Engel  _et al._ (1975), the nuclear part of which
can be written as a general energy density functional (EDF) of time-even
densities Dobaczewski and Dudek (1995):

${\cal
H}_{t}(r)=C_{t}^{\rho}\rho_{t}^{2}+C_{t}^{\rho\Delta\rho}\rho_{t}\Delta\rho_{t}+C_{t}^{\tau}\rho_{t}\tau_{t}+C_{t}^{J}\tensor{J}_{t}^{2}+C_{t}^{\rho\nabla
J}\rho_{t}\nabla\cdot\mathbf{J}_{t},$ (8)

where the subscript $t=(0,1)$ represents isoscalar and isovector densities,
respectively. The parameters of this EDF, $C_{t}^{\tau}$ for instance, model
the coupling between the particles and the nucleonic density in question (the
kinetic energy density, $\tau$, in this case). As it is usually done in modern
EDF optimization Chabanat  _et al._ (1997); Kortelainen  _et al._ (2010), we
can parametrize those couplings in terms of nuclear-matter properties (NMP)
plus the remaining coupling constants unconstrained by NMP:

$\left\\{\rho_{\text{c}},E^{\text{NM}}/A,K^{\text{NM}},a_{\text{sym}}^{\text{NM}},L_{\text{sym}}^{\text{NM}},M_{s}^{*},C_{t}^{\rho\Delta\rho},C_{t}^{\rho\nabla
J}\right\\}.$ (9)

This representation is primarily rooted in physical observables—like the
nuclear saturation density, $\rho_{c}$—and simplifies the selection of a
sensible range of values to explore in model calibration for DFT, and for
constructing the training bases for the RBM.

To test the RBM in extrapolation for DFT, we built a Lagrange RB of four
points spanning $L_{\text{sym}}^{\text{NM}}=[30,60]$ MeV and
$K^{\text{NM}}=[200,220]$ MeV while the other parameters remained at their
optimized UNEDF1 values Kortelainen  _et al._ (2012). The wave functions on
each shell (7 for neutrons and 6 for protons) were calculated using both the
exact solver and the RBM. Panel b) in Fig. 2 shows the performance of the RBM
when calculating the neutron skin thickness of 48Ca Hagen  _et al._ (2016), a
quantity particularly sensitive to the $L_{\text{sym}}^{\text{NM}}$ parameter.
The agreement between the emulated values and exact DFT results is excellent,
with an error of less than $0.8\%$ for all extrapolated values shown, even for
$L_{\text{sym}}^{\text{NM}}$ and $K^{\text{NM}}$ well outside the training
zone.

To test the limits of the RBM, the range of all ten available parameters in
Eq. (9) were widened well beyond what is reasonable for realistic NMP. We used
LHS to draw 50 training points to build a POD RB with $n=(2,4,8)$ and to
independently draw 500 testing points within the widened NMP ranges. As such,
several parameter combinations yielded convergence issues for the DFT solver,
but not for the RBM calculations, highlighting the capability of RBMs to
extrapolate into regions where exact solvers can experience numerical
instabilities. Even though the emulated results of non-converging test points
seemed reasonable, we consider their validation to be beyond the scope of this
work.

As Table 1 shows, for the stable parameter sets, the RBM reproduces single
nucleon energies well. This is particularly striking for the RB with only two
elements, which gives an error of about $0.6\%$ despite all ten parameters
being varied in the test sample. In terms of speed-up when calculating the 500
testing points, the RB with $n=(2,4,8)$ were $6$, $4$, and $2$ times faster
than the exact solver, respectively. We expect greater speed-ups to be
achieved by using techniques developed for nonlinear problems, such as
Barrault  _et al._ (2004); Grepl  _et al._ (2007), and by integrating the RBM
directly into the routine that calculates the DFT solutions, in a strategy
similar to the offline/online procedure often seen in the RBM literature
Quarteroni  _et al._ (2015). This integration could enable large scale model
UQ studies for a wide range of EDFs, an endeavor which, up to now, seemed
inaccessible. Furthermore, the RBM approach could also reduce the penalty of
using higher-dimensional solvers for systematic studies and UQ, calculations
previously limited to spherical and cylindrical symmetries.

We hope our results help spark the interest of the nuclear theory community in
RBMs. For this purpose, we created and will continue to update an online
resource rbm  to illustrate many of the concepts we discussed. The adoption of
recent developments on the choice of ansatz subspaces Benner  _et al._ (2015);
Ohlberger and Rave (2013); Nonino  _et al._ (2019), on error bounds and
convergence properties of RBMs Veroy  _et al._ (2003); Buffa  _et al._ (2012);
Cohen  _et al._ (2020), and on the computational efficiency for non-affine and
nonlinear problems Barrault  _et al._ (2004); Grepl  _et al._ (2007), to name
a few, could become key in reaching the full extent of what RBMs can offer.
Given the simplicity and flexibility of the Galerkin projection, and the PCA
diagnostic we showcase to test for low-dimensional manifolds, we believe that
RBMs have the potential to become standard tools for the emulation of
challenging problems in many-body nuclear physics.

##  I Acknowledgments

###### Acknowledgements.

While preparing this manuscript, we learned in Melendez (2021a) that Jordan
Melendez and collaborators were independently working on an overview of model
order reduction for nuclear systems Melendez (2021b); Melendez  _et al._
(2022). We are very grateful to Witek Nazarewicz and Frederi Viens for their
critical observations during the elaboration of this manuscript. We also thank
Ana Posada, Jorge Piekarewicz, and Daniel Phillips for their careful read of
the manuscript. This work was supported by the National Science Foundation
CSSI program under award number 2004601 (BAND collaboration), the U.S.
Department of Energy under Award Number DOE-DE-NA0003885 (NNSA, the
Stewardship Science Academic Alliances program), U.S. Department of Energy
(DE-SC0013365 and DE-SC0021152) and the Nuclear Computational Low-Energy
Initiative (NUCLEI) SciDAC-4 project (DE-SC0018083).

## References

  * Phillips  _et al._ (2021) D R Phillips, R J Furnstahl, U Heinz, T Maiti, W Nazarewicz, F M Nunes, M Plumlee, M T Pratola, S Pratt, F G Viens, and S M Wild, “Get on the BAND wagon: a Bayesian framework for quantifying model uncertainties in nuclear dynamics,” J. Phys. G Nucl. Part. Phys. 48, 072001 (2021).
  * Ekström  _et al._ (2019) Andreas Ekström, Christian Forssén, Christos Dimitrakakis, D Dubhashi, Håkan T Johansson, Azam Sheikh Muhammad, Hans Salomonsson, and Alexander Schliep, “Bayesian optimization in _ab initio_ nuclear physics,” J. Phys. G Nucl. Part. Phys. 46, 095101 (2019).
  * Melendez  _et al._ (2021a) JA Melendez, RJ Furnstahl, HW Grießhammer, JA McGovern, DR Phillips, and MT Pratola, “Designing optimal experiments: an application to proton Compton scattering,” Eur. Phys. J. A 57, 1–24 (2021a).
  * Giuliani and Piekarewicz (2021) PG Giuliani and J Piekarewicz, “From noise to information: The transfer function formalism for uncertainty quantification in reconstructing the nuclear density,” Phys. Rev. C 104, 024301 (2021).
  * Pratt  _et al._ (2015) Scott Pratt, Evan Sangaline, Paul Sorensen, and Hui Wang, “Constraining the Equation of State of Superhadronic Matter from Heavy-Ion Collisions,” Phys. Rev. Lett. 114, 202301 (2015).
  * McDonnell  _et al._ (2015) J. D. McDonnell, N. Schunck, D. Higdon, J. Sarich, S. M. Wild, and W. Nazarewicz, “Uncertainty Quantification for Nuclear Density Functional Theory and Information Content of New Measurements,” Phys. Rev. Lett. 114, 122501 (2015).
  * Wesolowski  _et al._ (2019) S Wesolowski, RJ Furnstahl, JA Melendez, and DR Phillips, “Exploring Bayesian parameter estimation for chiral effective field theory using nucleon–nucleon phase shifts,” J. Phys. G Nucl. Part. Phys. 46, 045102 (2019).
  * Lasseri  _et al._ (2020) Raphaël-David Lasseri, David Regnier, Jean-Paul Ebran, and Antonin Penon, “Taming Nuclear Complexity with a Committee of Multilayer Neural Networks,” Phys. Rev. Lett. 124, 162502 (2020).
  * Santner  _et al._ (2003) Thomas J Santner, Brian J Williams, William I Notz, and Brain J Williams,  _The Design and Analysis of Computer Experiments_, Vol. 1 (Springer, 2003).
  * Frame  _et al._ (2018) Dillon Frame, Rongzheng He, Ilse Ipsen, Daniel Lee, Dean Lee, and Ermal Rrapaj, “Eigenvector Continuation with Subspace Learning,” Phys. Rev. Lett. 121, 032501 (2018).
  * König  _et al._ (2020) S. König, A. Ekström, K. Hebeler, D. Lee, and A. Schwenk, “Eigenvector continuation as an efficient and accurate emulator for uncertainty quantification,” Phys. Lett. B 810, 135814 (2020).
  * Furnstahl  _et al._ (2020) R.J. Furnstahl, A.J. Garcia, P.J. Millican, and Xilin Zhang, “Efficient emulators for scattering using eigenvector continuation,” Phys. Lett. B 809, 135719 (2020).
  * Melendez  _et al._ (2021b) J.A. Melendez, C. Drischler, A.J. Garcia, R.J. Furnstahl, and Xilin Zhang, “Fast & accurate emulation of two-body scattering observables without wave functions,” Phys. Lett. B 821, 136608 (2021b).
  * Drischler  _et al._ (2021) C. Drischler, M. Quinonez, P.G. Giuliani, A.E. Lovell, and F.M. Nunes, “Toward emulating nuclear reactions using eigenvector continuation,” Phys. Lett. B 823, 136777 (2021).
  * Cohen-Tannoudji  _et al._ (1986a) Claude Cohen-Tannoudji, Bernard Diu, and Frank Laloe,  _Quantum Mechanics_ , Vol. 2 (John Wiley & Sons, 1986).
  * Kohn (1948) W. Kohn, “Variational Methods in Nuclear Collision Problems,” Phys. Rev. 74, 1763–1772 (1948).
  * Almroth  _et al._ (1978) B. O. Almroth, Perry Stern, and Frank A Brogan, “Automatic choice of global shape functions in structural analysis,” AIAA J. 16, 525–528 (1978).
  * Quarteroni  _et al._ (2015) Alfio Quarteroni, Andrea Manzoni, and Federico Negri,  _Reduced Basis Methods for Partial Differential Equations: An Introduction_, Vol. 92 (Springer, 2015).
  * Hesthaven  _et al._ (2016) Jan S Hesthaven, Gianluigi Rozza, Benjamin Stamm, _et al._ ,  _Certified Reduced Basis Methods for Parametrized Partial Differential Equations_, Vol. 590 (Springer, 2016).
  * Quarteroni  _et al._ (2011) Alfio Quarteroni, Gianluigi Rozza, and Andrea Manzoni, “Certified reduced basis approximation for parametrized partial differential equations and applications,” J. Math. Ind. 1, 1–49 (2011).
  * Nguyen  _et al._ (2010) N. C. Nguyen, G. Rozza, D. B. P. Huynh, and A. T. Patera, “Reduced Basis Approximation and a Posteriori Error Estimation for Parametrized Parabolic PDEs: Application to Real-Time Bayesian Parameter Estimation,” in  _Large‐Scale Inverse Problems and Quantification of Uncertainty_ (John Wiley & Sons, Ltd, 2010) Chap. 8, pp. 151–177.
  * Field  _et al._ (2011) Scott E Field, Chad R Galley, Frank Herrmann, Jan S Hesthaven, Evan Ochsner, and Manuel Tiglio, “Reduced Basis Catalogs for Gravitational Wave Templates,” Phys. Rev. Lett. 106, 221102 (2011).
  * Milani  _et al._ (2008) Roberto Milani, Alfio Quarteroni, and Gianluigi Rozza, “Reduced basis method for linear elasticity problems with many parameters,” Comput. Methods Appl. Mech. Eng. 197, 4812–4829 (2008).
  * Thompson and Nunes (2009) Ian J Thompson and Filomena M Nunes,  _Nuclear Reactions for Astrophysics: Principles, Calculation and Applications of Low-Energy Reactions_ (Cambridge University Press, 2009).
  * Jolliffe (2002) Ian T Jolliffe,  _Principal Component Analysis_ (Springer, 2002).
  * Blum  _et al._ (2020) Avrim Blum, John Hopcroft, and Ravindran Kannan,  _Foundations of Data Science_ (Cambridge University Press, 2020).
  * Nonino  _et al._ (2019) Monica Nonino, Francesco Ballarin, Gianluigi Rozza, and Yvon Maday, “Overcoming slowly decaying kolmogorov n-width by transport maps: application to model order reduction of fluid dynamics and fluid–structure interaction problems,”  (2019), arXiv:1911.06598 [math.NA] .
  * Rawitscher  _et al._ (2018) George Rawitscher, Victo dos Santos Filho, and Thiago Carvalho Peixoto, “Galerkin and Collocation Methods,” in  _An Introductory Guide to Computational Methods for the Solution of Physics Problems_ (Springer, 2018) pp. 17–31.
  * Fletcher (1984) Clive AJ Fletcher,  _Computational Galerkin Methods_ (Springer, 1984).
  * Buffa  _et al._ (2012) Annalisa Buffa, Yvon Maday, Anthony T Patera, Christophe Prud’homme, and Gabriel Turinici, “ _A priori_ convergence of the Greedy algorithm for the parametrized reduced basis method,” ESAIM: Math. Model. Numer. Anal. 46, 595–603 (2012).
  * Prud’homme  _et al._ (2001) C. Prud’homme, D. V. Rovas, K. Veroy, L. Machiels, Y. Maday, A. T. Patera, and G. Turinici, “Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods ,” J. Fluids Eng. 124, 70–80 (2001).
  * Lucchese (1989) Robert R. Lucchese, “Anomalous singularities in the complex Kohn variational principle of quantum scattering theory,” Phys. Rev. A 40, 6879–6885 (1989).
  * Singh and Stauffer (1974) SR Singh and AD Stauffer, “A unified formulation of the variational methods in scattering theory,” Il Nuovo Cimento B (1971-1996) 22, 139–152 (1974).
  * Cohen-Tannoudji  _et al._ (1986b) Claude Cohen-Tannoudji, Bernard Diu, and Frank Laloe,  _Quantum Mechanics_ , Vol. 1 (John Wiley & Sons, 1986).
  * Thompson  _et al._ (1977) DR Thompson, M LeMere, and YC Tang, “Systematic investigation of scattering problems with the resonating-group method,” Nucl. Phys. A 286, 53–66 (1977).
  * Gross (1961) E. P. Gross, “Structure of a quantized vortex in boson systems,” Il Nuovo Cimento 20, 454–477 (1961).
  * Pitaevskii (1961) L. P. Pitaevskii, “Vortex Lines in an Imperfect Bose Gas,” JETP 13, 451 (1961).
  * Pichi  _et al._ (2020) Federico Pichi, Annalisa Quaini, and Gianluigi Rozza, “A Reduced Order Modeling Technique to Study Bifurcating Phenomena: Application to the Gross–Pitaevskii Equation,” SIAM J. Sci. Comput. 42, B1115–B1135 (2020).
  * Carr  _et al._ (2000a) L. D. Carr, Charles W. Clark, and W. P. Reinhardt, “Stationary solutions of the one-dimensional nonlinear Schrödinger equation. I. Case of repulsive nonlinearity,” Phys. Rev. A 62, 063610 (2000a).
  * Carr  _et al._ (2000b) L. D. Carr, Charles W. Clark, and W. P. Reinhardt, “Stationary solutions of the one-dimensional nonlinear Schrödinger equation. II. Case of attractive nonlinearity,” Phys. Rev. A 62, 063611 (2000b).
  * Carr  _et al._ (2001) L. D. Carr, C. W. Clark, and W. P. Reinhardt, “Expected and Unexpected Solutions to the Stationary One-Dimension Nonlinear Schödinger Equation,” Int. J. Mod. Phys. B 15, 1663 (2001).
  * Torres-Vega (2017) G Torres-Vega, “Some real, periodic stationary solutions of the one-dimensional nonlinear Schrödinger equation for constant potentials,” J. Phys. Conf. Ser. 792, 012054 (2017).
  * Mckay  _et al._ (2000) M. D. Mckay, R. J. Beckman, and W. J. Conover, “A comparison of three methods for selecting values of input variables in the analysis of output from a computer code,” Technometrics 42, 55–61 (2000).
  * Veroy  _et al._ (2003) Karen Veroy, Christophe Prud’Homme, Dimitrios Rovas, and Anthony Patera, “A Posteriori Error Bounds for Reduced-Basis Approximation of Parametrized Noncoercive and Nonlinear Elliptic Partial Differential Equations,” in  _16th AIAA Comput. Fluid Dyn. Conf._ (2003) p. 3847.
  * Haasdonk and Ohlberger (2008) Bernard Haasdonk and Mario Ohlberger, “Reduced basis method for finite volume approximations of parametrized linear evolution equations,” ESAIM: M2AN 42, 277–302 (2008).
  * Sarkar and Lee (2021) Avik Sarkar and Dean Lee, “Self-learning Emulators and Eigenvector Continuation,” (2021), arXiv:2107.13449 [nucl-th] .
  * Schunck (2019) Nicolas Schunck, ed.,  _Energy Density Functional Methods for Atomic Nuclei_, 2053-2563 (IOP Publishing, 2019).
  * Cancès  _et al._ (2007) Éric Cancès, Claude Le Bris, Yvon Maday, Ngoc Nguyen, Anthony Patera, and George Pau, “Feasibility and competitiveness of a reduced basis approach for rapid electronic structure calculations in quantum chemistry,” in  _CRM Proceedings and Lecture Notes_ (American Mathematical Society, 2007) pp. 15–47.
  * Lin  _et al._ (2012) Lin Lin, Jianfeng Lu, Lexing Ying, and E Weinan, “Adaptive local basis set for Kohn–Sham density functional theory in a discontinuous Galerkin framework I: Total energy calculation,” J. Comput. Phys. 231, 2140–2154 (2012).
  * Zhang  _et al._ (2017) Gaigong Zhang, Lin Lin, Wei Hu, Chao Yang, and John E Pask, “Adaptive local basis set for Kohn–Sham density functional theory in a discontinuous Galerkin framework II: Force, vibration, and molecular dynamics calculations,” J. Comput. Phys. 335, 426–443 (2017).
  * Skyrme (1956) T. H. R. Skyrme, “CVII. The nuclear surface,” Philos. Mag. 1, 1043–1054 (1956).
  * Skyrme (1958) T.H.R. Skyrme, “The effective nuclear potential,” Nucl. Phys. 9, 615–634 (1958).
  * Vautherin and Brink (1972) D. Vautherin and D. M. Brink, “Hartree-Fock Calculations with Skyrme’s Interaction. I. Spherical Nuclei,” Phys. Rev. C 5, 626–647 (1972).
  * Engel  _et al._ (1975) Y.M. Engel, D.M. Brink, K. Goeke, S.J. Krieger, and D. Vautherin, “Time-dependent Hartree-Fock theory with Skyrme’s interaction,” Nucl. Phys. A 249, 215–238 (1975).
  * Dobaczewski and Dudek (1995) J. Dobaczewski and J. Dudek, “Time-odd components in the mean field of rotating superdeformed nuclei,” Phys. Rev. C 52, 1827–1839 (1995).
  * Chabanat  _et al._ (1997) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, “A Skyrme parametrization from subnuclear to neutron star densities,” Nucl. Phys. A 627, 710–746 (1997).
  * Kortelainen  _et al._ (2010) M. Kortelainen, T. Lesinski, J. Moré, W. Nazarewicz, J. Sarich, N. Schunck, M. V. Stoitsov, and S. Wild, “Nuclear energy density optimization,” Phys. Rev. C 82, 024313 (2010).
  * Kortelainen  _et al._ (2012) M. Kortelainen, J. McDonnell, W. Nazarewicz, P.-G. Reinhard, J. Sarich, N. Schunck, M. V. Stoitsov, and S. M. Wild, “Nuclear energy density optimization: Large deformations,” Phys. Rev. C 85, 024304 (2012).
  * Hagen  _et al._ (2016) G Hagen, A Ekström, C Forssén, GR Jansen, W Nazarewicz, T Papenbrock, KA Wendt, Sonia Bacca, N Barnea, B Carlsson, _et al._ , “Neutron and weak-charge distributions of the 48Ca nucleus,” Nat. Phys. 12, 186–190 (2016).
  * Barrault  _et al._ (2004) Maxime Barrault, Yvon Maday, Ngoc Cuong Nguyen, and Anthony T Patera, “An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations,” C. R. Math. 339, 667–672 (2004).
  * Grepl  _et al._ (2007) Martin A Grepl, Yvon Maday, Ngoc C Nguyen, and Anthony T Patera, “Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations,” ESAIM: Math. Model. Numer. Anal. 41, 575–605 (2007).
  * (62) https://kylegodbey.github.io/nuclear-rbm.
  * Benner  _et al._ (2015) Peter Benner, Serkan Gugercin, and Karen Willcox, “A survey of projection-based model reduction methods for parametric dynamical systems,” SIAM Review 57, 483–531 (2015).
  * Ohlberger and Rave (2013) Mario Ohlberger and Stephan Rave, “Nonlinear reduced basis approximation of parameterized evolution equations via the method of freezing,” C. R. Math. 351, 901–906 (2013).
  * Cohen  _et al._ (2020) Albert Cohen, Wolfgang Dahmen, Ronald DeVore, and James Nichols, “Reduced Basis Greedy Selection Using Random Training Sets,” ESAIM: M2AN 54, 1509–1524 (2020).
  * Melendez (2021a) Jordan Melendez, “Variational principles & efficient subspace emulation. ISNET Conference December 2021,”  (2021a).
  * Melendez (2021b) Jordan Melendez, “Variational principles & efficient subspace emulation,”  (2021b).
  * Melendez  _et al._ (2022) J. Melendez, C. Drischler, Xilin Zhang, A. Garcia, and RJ Furnstahl, “In preparation,” (2022).
  * Ascher  _et al._ (1995) Uri M Ascher, Robert MM Mattheij, and Robert D Russell,  _Numerical solution of boundary value problems for ordinary differential equations_ (SIAM, 1995).

##  Appendix A Supplementary Material

##  Appendix B Equivalence of the RBM and Eigenvector Continuation under the
Kohn Variational Principle

For the sake of conciseness, we only treat the case of eigenvector
continuation (EC) applied to single-channel scattering, as done in Furnstahl
_et al._ (2020). Within this context, we aim at finding approximate solutions
to a differential equation of the form Thompson and Nunes (2009):

$F_{\alpha}(\phi)=\left(-\frac{d^{2}}{dr^{2}}+\frac{\ell(\ell+1)}{r^{2}}+U(r,\alpha)-p^{2}\right)\phi(r)=0,$
(S1)

where $U(r,\alpha)$ represent the potential the particle interacts with, which
we assume changes smoothly with its parameter set $\alpha$, $\ell$ is the
angular momentum quantum number and $p$ relates to the asymptotic momentum.
Let us assume that the solution to this equation is subject to the boundary
conditions $\phi(r=0)=0$, and

$\phi(r)\underset{r\rightarrow\infty}{\longrightarrow}\frac{1}{p}\sin\left(pr-\ell\frac{\pi}{2}\right)+\tau\cos\left(pr-\ell\frac{\pi}{2}\right).$
(S2)

Note that Eq. (S2) imposes a normalization condition on $\phi$: the
coefficient accompanying the sine function must equal $1/p$.

A straightforward application of the reduced basis method (RBM) as discussed
in the main document leads to the choice of an approximate function:

$\hat{\phi}_{\alpha}=\sum_{k=1}^{n}a_{k}\phi_{k},\quad\text{with}\quad\sum_{k=1}^{n}{a_{k}}=1,$
(S3)

where the $\phi_{k}$ are solutions to $F_{\alpha_{k}}\phi=0$ with the correct
boundary conditions. Following the discussion on the main text, we eliminate
the redundancy of the coefficients $a_{k}$ created by the boundary conditions
by explicitly writing one of the them in terms of the others. Without loss of
generality, we let $a_{1}=1-\sum_{k=2}^{n}a_{k}$, obtaining:

$\hat{\phi}_{\alpha}=\phi_{1}+\sum_{k=2}^{n}a_{k}(\phi_{k}-\phi_{1}).$ (S4)

With this rewriting, it becomes explicit that we only need to determine
$(n-1)$ coefficients to obtain the approximate solution, which is part of the
affine space spanned by $(\phi_{k}-\phi_{1})$. Following a similar reasoning
to the main text, we choose $\psi_{k}=({\phi}_{k}-\phi_{1})$ for
$k=2,\dots,n$, resulting in the following equations for the remaining $a_{k}$:

$\displaystyle\begin{split}&\langle\psi_{k}|F_{\alpha}(\hat{\phi}_{\alpha})\rangle=\langle{\phi}_{k}-\phi_{1}|F_{\alpha}(\phi_{1})\rangle\
+\\\
&\sum_{j=2}^{n}a_{j}\langle{\phi}_{k}-\phi_{1}|F_{\alpha}({\phi}_{j}-\phi_{1})\rangle=0,\end{split}$
(S5)

for $k=2\dots,n,$ where the inner product is defined as
$\langle\psi|\phi\rangle=\int_{0}^{\infty}\psi(r)\phi(r)dr$ (without the usual
complex conjugation).

We can rewrite Eq. (S5) in a more compact form by reintroducing $a_{1}$, and
the normalization condition. This results in the equation set:

$\displaystyle\sum_{j=1}^{n}a_{j}\langle{\phi}_{k}-\phi_{1}|F_{\alpha}(\phi_{j})\rangle=0,\quad\text{for}\quad
k=2,\dots,n;$ (S6) $\displaystyle\text{and}\quad a_{1}=1-\sum_{j=2}^{n}a_{j}.$
(S7)

###  B.1 Approximation through the Kohn Variational Principle

The ($K$-matrix) Kohn variational principle (KVP) states that the solution to
Eq. (S1) with the asymptotic behavior (S2), is a stationary point for the
functional:

$\beta[\phi]=\tau[\phi]-\int_{0}^{\infty}dr\phi(r)F_{\alpha}\left(\phi(r)\right),$
(S8)

where $\tau[\phi]$ extracts its value from the asymptotic behavior of $\phi$,
that is, the cosine coefficient in Eq. (S2).

The Kohn variational method is detailed in the supplemental material of Ref.
Furnstahl  _et al._ (2020). It utilizes a trial function $\hat{\phi}_{\alpha}$
constructed exactly as in Eq. (S3), and it finds a stationary point of the
functional (S8) in terms of the coefficients $a_{k}$. After using the method
of Lagrange multipliers to enforce the normalization condition, it is found
that the following equation set describes the stationary point:

$\tau_{k}-\lambda-\sum_{j=1}^{n}(\Delta U_{jk}+\Delta
U_{kj})a_{j}=0\quad\text{for}\quad k=1\dots,n;$ (S9)

where $\tau_{k}$ is the cosine coefficient in Eq. (S2) associated with each
$\phi_{k}$, and the matrix $\Delta U_{jk}$ is a shorthand notation for the
inner products:

$\Delta U_{jk}=\langle\phi_{j}|F_{\alpha}(\phi_{k})\rangle.$ (S10)

These equations, together with the normalization condition
$\sum_{k=1}^{n}{a_{k}}=1$, can be used to find the $n$ coefficients $a_{k}$,
plus the Lagrange multiplier $\lambda$.

###  B.2 Proof of equivalence between the methods

To compare Eq. (S9) to Eq. (S6), we need to rewrite Eq. (S9) in terms of
$\Delta U_{kj}$ only, eliminating both the $\tau_{k}$ and the Lagrange
multiplier $\lambda$. We can relate the elements of $\Delta U_{jk}$ to their
transposes $\Delta U_{kj}$ by integration by parts. Note that

$-\int_{0}^{\infty}\phi_{j}\phi^{\prime\prime}_{k}dr=-\left.\left(\phi_{j}\phi_{k}^{\prime}-\phi_{j}^{\prime}\phi_{k}\right)\right|_{0}^{\infty}-\int_{0}^{\infty}\phi_{j}^{\prime\prime}\phi_{k}dr,$
(S11)

where the boundary term can be evaluated through the boundary conditions (S2)
to be $\tau_{k}-\tau_{j}$. Therefore, via integration by parts, we can make
$F_{\alpha}$ act on $\phi_{j}$ in Eq. (S10) to obtain:

$\Delta U_{jk}=\Delta U_{kj}+(\tau_{k}-\tau_{j}).$ (S12)

Using this result we can eliminate $\Delta U_{jk}$ from Eq. (S9) and obtain:

$\displaystyle\begin{split}\tau_{k}-\lambda-\sum_{j=1}^{n}(2\Delta
U_{kj})a_{j}-\sum_{j=1}^{n}a_{j}(\tau_{k}-\tau_{j})=0\\\
-\lambda-\sum_{j=1}^{n}(2\Delta
U_{kj})a_{j}+\sum_{j=1}^{n}a_{j}\tau_{j}=0,\end{split}$ (S13)

for $k=1,\dots,n$; where we used the fact that the $a_{k}$ sum to unity to
cancel the $\tau_{k}$ term.

Next, we can eliminate the Lagrange multiplier $\lambda$ and the sum of
$a_{j}\tau_{j}$ by subtracting equations for different $k$. Without loss of
generality, we can subtract all equations to the equation corresponding to
$k=1$, resulting in:

$\sum_{j=1}^{n}(\Delta U_{kj}-\Delta U_{1j})a_{j}=0\quad\text{for}\quad
k=2\dots,n.$ (S14)

Finally, using the definition of $\Delta
U_{kj}=\langle\phi_{k}|F_{\alpha}(\phi_{j})\rangle$ shows the equivalence with
Eq. (S6), completing the proof.

It is interesting to note that the second term in Eq. (S2) can be thought as a
first-order correction to the emulated $\tau$. This correction adds a factor
proportional to the residual $F_{\alpha}(\hat{\phi}_{\alpha})$, in a similar
spirit as Newton’s method Ascher  _et al._ (1995). The ‘proportionally factor’
(or transfer function Giuliani and Piekarewicz (2021)) connecting the residual
$|F_{\alpha}(\hat{\phi})\rangle$ and $\tau$ is in this case the negative of
the exact solution $-\langle\phi|$. This in turn, can be approximated by the
emulated solution $-\langle\hat{\phi}|$, obtaining the correction:

$-\langle\hat{\phi}|F_{\alpha}(\hat{\phi})\rangle=-\int_{0}^{\infty}dr\phi(r)F_{\alpha}\left(\phi(r)\right).$
(S15)

Such first order correction, together with the POD+Greedy algorithm we used in
the main text, highlight the fact that the residual $F_{\alpha}(\hat{\phi})$
is information rich and can be useful to improve the performance of the RBM
emulator.

###  B.3 Extension to the generalized KVP

This proof can be extended to show the equivalence between the RBM and the
generalized KVP, as expressed in Drischler  _et al._ (2021); Lucchese (1989).
The generalized KVP applies to different normalizations of the asymptotic
behavior of candidate solutions for Eq. (S1):

$\phi(r)\underset{r\rightarrow\infty}{\longrightarrow}\frac{1}{\mathcal{N}}\phi^{(0)}(r)+\frac{L}{\mathcal{N}}\phi^{(1)}(r),$
(S16)

where $\mathcal{N}$ is an arbitrary scaling constant, and $\phi^{(0)}$ and
$\phi^{(1)}$ are two linearly independent asymptotic solutions, which can be
written as:

$\begin{bmatrix}\phi^{(0)}(r)\\\
\phi^{(1)}(r)\end{bmatrix}=\begin{bmatrix}u_{00}&u_{01}\\\
u_{10}&u_{11}\end{bmatrix}\begin{bmatrix}\sin\left(pr-\ell\frac{\pi}{2}\right)\\\
\cos\left(pr-\ell\frac{\pi}{2}\right)\end{bmatrix},$ (S17)

for a non-singular (generally complex) matrix $\\{\mathbf{u}\\}_{ij}=u_{ij}$.

The Generalized KVP states that the functional

$\beta_{\mathbf{u}}[\phi]=\frac{1}{\mathcal{N}}L[\phi]-\frac{\mathcal{N}}{p\det{\mathbf{u}}}\int_{0}^{\infty}dr\phi(r)F_{\alpha}\left(\phi(r)\right),$
(S18)

is stationary in the solutions to equation (S1) with the general asymptotic
behavior (S16). Proceeding in a manner similar to the $K-$matrix KVP case, an
approximation of the form (S3) can be found using the method of Lagrange
multipliers with the constraint $\sum_{k=1}^{n}a_{k}=1$. This results in the
equations:

$\frac{L_{k}}{\mathcal{N}}-\lambda-\frac{\mathcal{N}}{p\det{\mathbf{u}}}\sum_{j=1}^{n}(\Delta
U_{jk}+\Delta U_{kj})a_{j}=0,$ (S19)

for $k=1\dots,n$. Once again, $\Delta
U_{jk}=\langle\phi_{j}|F_{\alpha}(\phi)_{k}\rangle=\int_{0}^{\infty}dr\phi_{j}(r)F_{\alpha}\left(\phi_{k}(r)\right)$.
Note that in this definition for the product there is no complex conjugate on
either entry.

An integration by parts argument similar to Eq. (S11), leads to the relation:

$\Delta U_{jk}=\Delta
U_{kj}+\frac{p\det{\mathbf{u}}}{\mathcal{N}^{2}}(L_{k}-L_{j}).$ (S20)

The rest of the proof is analogous to Eqs. (S13-S14), where we substitute the
transposed array, and eliminate the Lagrange multiplier to obtain equation
(S6).

It is possible that a different normalization condition together with a
different choice of the reduced basis and projecting functions, for example a
POD basis informed on $N>n$ training points, would be more robust against the
so-called Kohn anomalies Kohn (1948). Such investigation is left for a future
work.

##  Appendix C RBMs on a set of coupled equations

As we described in the main text, the RBM can be applied to the case with a
set of $m$ coupled equations:

$F_{\alpha}^{(i)}(\phi^{(1)}_{\alpha},\dots,\phi^{(m)}_{\alpha})=0,\quad\text{for}\quad
i=1,\dots,m;$ (S21)

by approximating each $\phi_{\alpha}^{(i)}$ as a linear combination of their
corresponding solutions for different values of $\alpha$:

$\displaystyle\begin{split}&\hat{\phi}^{(i)}_{\alpha}=\phi_{0}^{(i)}+\sum_{k=1}^{n}a^{(i)}_{k}\phi_{k}^{(i)},\\\
&\text{for}\quad k=1,\dots,n;\quad i=1,\dots,m.\end{split}$ (S22)

By selecting $\psi_{k}^{(i)}$ as the generators of the affine spaces for each
$\hat{\phi}^{(i)}_{\alpha}$, we obtain the $n\times m$ Galerkin equations:

$\langle\psi_{k}^{(i)}|F^{(i)}_{\alpha}(\phi^{(1)}_{\alpha},\dots,\phi^{(m)}_{\alpha})\rangle=0,$
(S23)

for $k=1,\dots,n;\ \text{and}\ i=1,\dots,m;$ to obtain the coefficients
$a_{k}^{(i)}$.

In the case of coupled eigenvalue-eigenvector systems, we can proceed as in
the case of a single equation, by substituting the $m$ eigenvalues
$\lambda^{(i)}_{\alpha}$ for approximate values
$\hat{\lambda}^{(i)}_{\alpha}$, and enforcing $m$ normalization conditions, in
accordance with the requirements of the problem at hand.

##  Appendix D Details about the numerical results

###  D.1 Decay of singular values

In this section, we give details on the construction of the singular values
$\sigma_{k}$ from the singular value decomposition (SVD) showed in Fig. 1 of
the main text. Four problems were considered: the infinite well (IW), the
single channel 2-body scattering with a Minnesota potential at fixed energy
(SM$E$) and varying energy (SM$E^{*}$), the Gross-Pitaevskii equation (GP),
and 48Ca under density functional theory (DFT). Fig S1 shows, for all problems
considered, solutions for 10 different values of their parameters (left
column), and the first 4 principal components out of a sample of 40 parameters
in each case (right column).

The IW problem consists of the 1D quantum Hamiltonian of a particle trapped in
an infinite well (IW) Cohen-Tannoudji  _et al._ (1986b) where $\alpha$
controls the location of the well:

$V(x,\alpha)=\begin{cases}0,&\alpha<x<\alpha+1,\\\
\infty,&\text{otherwise.}\end{cases}$ (S24)

The ground state solutions to this Hamiltonian are wave functions of the form
$\phi_{\alpha}=\sqrt{2}\sin{[\pi(x-\alpha)]}$ for $\alpha<x<\alpha+1$ and zero
otherwise. The singular values $\sigma_{k}$ for the IW showed in Fig. 1 of the
main text were obtained by sampling 40 values of $\alpha$ in the range
$[-5,5]$ using Latin hypercube sampling (LHS), and performing SVD on the set
of 40 solutions. The $\sigma_{k}$ do not decay exponentially, and as can be
seen in Fig. S1 the first 4 principal components are unable to capture the
variability of the set of solutions. As mentioned in the main text, extensions
to the basic RBM Nonino  _et al._ (2019), such as allowing the training
functions to be re-scaled and shifted in their $x$ domain, are able to tackle
these issues. In this particular example, taking advantage of the symmetry of
the problem, i.e., $\phi_{\alpha}(x+\alpha)=\phi_{\alpha=0}(x)$, would lead to
all $\sigma_{k}$ to be zero for $k\geq 2$.

Both SM$E$ and SM$E^{*}$ problems consist of the single-channel ${}^{1}S_{0}$
nucleon-nucleon scattering Hamiltonian Thompson and Nunes (2009). The
Minnesota potential Thompson  _et al._ (1977) with parameters
$[V_{0R},V_{0S}]$ is used for the interaction:

$V(r,V_{0R},V_{0S})=V_{0R}e^{-1.487r^{2}}+V_{0S}e^{-0.465r^{2}},$ (S25)

where the other two non-linear parameters were fixed. In both cases, we make
the change of variables $s=pr$ and the scattering Hamiltonian takes the form:

$\left(-\frac{d^{2}}{ds^{2}}+\frac{\ell(\ell+1)}{s^{2}}+\tilde{U}(s,\alpha,p)-1\right)\phi_{\alpha}(s)=0,$
(S26)

where the potential $U(r,\alpha)=V(r,\alpha)2\mu$ is now momentum dependent:
$\tilde{U}(s,\alpha,p)=U(s/p,\alpha)/p^{2}$. In both cases, (SM$E$ and
SM$E^{*}$) 40 parameters where obtained by a LHS in the range
$V_{0R}=[100,300]$ MeV and $V_{0S}=[-200,0]$ MeV, following Furnstahl  _et
al._ (2020). The singular values showed in Fig. 1 of the main text were
obtained by performing SVD on the set of 40 solutions. In the case of SM$E$,
all 40 solutions shared the same energy in the center of mass $E=50$ MeV,
while for SM$E^{*}$ the energies where equispaced in the range $E=[20,80]$
MeV.

Figure S1: Examples of solutions (left column) and principal components (right
column) for all the problems discussed in Fig. 1 in the main text: panels a)
and b) correspond to the Infinite Well; panels c) and d) to the Scattering
with the Minnesota potential at a fixed energy; panels e) and f) to the
Scattering with the Minnesota potential with a variable energy; panels g) and
h) to the Gross-Pitaevskii equation; and panels i) and j) to 48Ca under DFT.
For panels i) and j) we chose the ${}^{2}S_{1/2}$ neutron wave function as a
showcase, out of the total of 13 distinct levels for protons and neutrons. The
left column shows exact solutions for 10 different values of the respective
sampled parameters, while the right column shows the first 4 principal
components obtained from a SVD analysis on a set of 40 exact solutions.

The GP and DFT cases are explained in detail in the main text. The ranges for
the parameters in both cases correspond to the ones used in Table I of the
main text. For GP, a set of 40 values of the parameters $[q,\kappa]$ were
obtained by LHS in the range $q\in[0,30]$ and $\kappa\in[5,30]$. For DFT, 50
values of the parameters were obtained with a LHS across the parameter ranges
shown in Table S1.

###  D.2 1-D Gross-Pitaevskii equation with a harmonic trapping potential

The four training points $\phi_{k}$ used in the Lagrange basis for the results
of Fig. 2 a) in the main text are :
$[q,\kappa]=\\{[0,1],[0,5],[0.5,1],[0.5,5]\\}$. Fig. S2 shows the training
points for the Lagrange and POD RB, as well as the 500 testing points used for
the results shown in Table I in the main text. Fig. S3 shows the construction
on the POD+Greedy basis also used for the results shown in Table I in the main
text. Fig. S4 shows the 20 exact solutions $\phi_{q,\kappa}(x)$ selected by
the Greedy algorithm, as well as the first four principal components of this
set.

Define starting parameters $\alpha_{1}$, $n_{2}$, $N$, $N_{2}$

Find $\tilde{\phi}_{1}$ s.t. $F_{\alpha_{1}}(\tilde{\phi}_{1})=0$

$ES\leftarrow\\{\tilde{\phi}_{1}\\}$

for  _$i=2,\dots,N$_ do

$RB\leftarrow POD(ES,n_{2})$

Draw $N_{2}$ parameters with LHS:
$A=\\{\hat{\alpha}_{1},\dots,\hat{\alpha}_{N_{2}}\\}$

Use the RBM with the RB to find $\hat{\phi}_{\hat{\alpha}}$ for each element
in $A$

$\alpha_{i}\leftarrow\operatorname*{arg\,max}_{\hat{\alpha}\in
A}||F_{\hat{\alpha}}(\hat{\phi}_{\hat{\alpha}})||^{2}$

Find $\tilde{\phi}_{i}$ s.t. $F_{\alpha_{i}}(\tilde{\phi}_{i})=0$

$ES$ $\leftarrow ES\cup\\{\tilde{\phi}_{i}\\}$

end for

return $ES$

Algorithm 1 POD+Greedy scheme Table S1: Ranges for the 10 parameters used to generate the DFT results in Table I in the main text. | Min | Max | Units
---|---|---|---
$\rho_{c}$ | 0.14 | 0.18 | fm-3
E${}^{\text{NM}}$/A | -16.5 | -14.5 | MeV
K${}^{\text{NM}}$ | 160 | 260 | MeV
a${}^{\text{NM}}_{\text{sym}}$ | 26 | 32 | MeV
L${}^{\text{NM}}_{\text{sym}}$ | 20 | 180 | MeV
M${}^{*}_{s}$ | 0.7 | 1.4 |
$C_{0}^{\rho\Delta\rho}$ | -55 | -40 | MeV fm5
$C_{1}^{\rho\Delta\rho}$ | -165 | -90 | MeV fm5
$C_{0}^{\rho\nabla J}$ | -105 | -55 | MeV fm5
$C_{1}^{\rho\nabla J}$ | -50 | -15 | MeV fm5

####  D.2.1 Description of the POD+Greedy scheme used for GP

The POD+Greedy scheme used for the results of Table I of the main text
consists on iteratively constructing a set of $N$ exact solutions ($ES$) with
a (weak) Greedy algorithm Quarteroni  _et al._ (2015) informed by the
residuals of a POD RB of dimension $n_{2}$ derived from the set of exact
solutions at each step. To set up the algorithm, let
$POD(\\{\tilde{\phi}_{l}\\}_{l=1}^{N},n)$ be a function that returns a
normalized basis constructed with first $n$ principal components of the set of
solutions $\\{\tilde{\phi}_{l}\\}_{l=1}^{N}$ if $N\geq n$, and returns
$POD(\\{\tilde{\phi}_{l}\\}_{l=1}^{N},N)$ if $N<n$. This function can be used
to construct a POD basis of size up to $n$ with a set of solutions. The Greedy
strategy we used is summarized in Algorithm 1. The desired POD bases were
constructed by running $POD(RB,n)$ on the output of this algorithm with
$\alpha_{1}=[q_{1},\kappa_{1}]=[15,17.5]$, $n_{2}=10$, $N=20$, and
$N_{2}=100$.

###  D.3 Spherical Nuclear Density Functional Theory

The four training points, $\phi_{k}$, used in Fig. 2 b) of the main text only
varied the $K^{\text{NM}}$ and $L_{\text{sym}}^{\text{NM}}$ parameters, with
the rest taken to be the standard UNEDF1 optimal parameters Kortelainen  _et
al._ (2012). The four values of $K^{\text{NM}}$ and
$L_{\text{sym}}^{\text{NM}}$, in MeV, are:

$\displaystyle
K^{\text{NM}}=200,L_{\text{sym}}^{\text{NM}}=30;K^{\text{NM}}=200,L_{\text{sym}}^{\text{NM}}=60;$
$\displaystyle
K^{\text{NM}}=220,L_{\text{sym}}^{\text{NM}}=30;K^{\text{NM}}=220,L_{\text{sym}}^{\text{NM}}=60.$

Table S1 shows the parameter ranges used for the LHS for the DFT results in
Table I of the main text. Both the 500 testing points and the $N=20$ exact
evaluations used to build the POD RB were independently drawn by LHS on these
ranges.

Figure S2: Parameter values used for the results of the GP equation in Table I
of the main text, all drawn by Latin Hypercube Sampling Mckay  _et al._
(2000). Points in blue show the 500 testing samples. Points in red, orange,
and green show the parameters for the Lagrange RB built with 2, 4, and 8 exact
calculations, respectively. Points in black show the 20 parameters used to
build the three POD RBs with $n=(2,4,8)$. Figure S3: Construction of the
POD+Greedy RB for the results in Table I of the main text. The black stars
show the $N$ exact calculations made at each stage to build the POD+Greedy
basis. The four panels show the stages for $N=(1,3,6,20)$. The red points in
each panel are the 100 LHS draws on which the norm of the residual is
maximized each time to add a new exact calculation to the POD+Greedy RB. This
random sampling is repeated on each step (no two panels share the same red
points). Points in blue and green on each panel show the 10 locations where
the residual is maximized or minimized, representing the regions where the
emulator is performing poorly and adequately, respectively. The green box
shows the limits for the LHS: $q\in[0,30]$ and $\kappa\in[5,30]$. Figure S4:
Panel a) shows the 20 exact solutions $\phi_{q,\kappa}(x)$ obtained by the
Greedy algorithm in Fig. S3 and used to construct the POD+Greedy basis for the
results of Table I in the main text. Panel b) shows the first four principal
components of this set which constitute the RB used for the POD+Greedy results
with $n=4$ in Table I in the main text.