Discovering reduced order model equations of many-body quantum systems using genetic programming: a technical report

Abstract

In this technical report we present first results of applying Genetic Programming (GP) to obtain equations for constructing reduced order models of quantum systems, with a particular interest in nuclear Density Functional Theory. We employ the reduced basis method to obtain reduced coordinates as the amplitudes of the reduced basis, and use GP to avoid the need of constructing the reduced equations through, for example, a Galerkin projection. The reduced order models constructed through GP show excellent accuracy and speed performance, including extrapolations in the controlling parameters, and show promise as an effective method for emulating computationally demanding calculations in nuclear physics.

I Introduction

Computational models for many-body quantum systems have quickly grown in complexity as theories get refined, more experimental data becomes available, and computers become more powerful. In fields with practical applications, such as nuclear physics, the demand for more accurate descriptions of nature, together with a push for predictions with well-quantified uncertainties [1, 2], has driven various efforts in the field to create surrogate models that can learn low-dimensional representations of the underlying system equations [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. Several of the developed surrogate models are based on model reduction approaches [30, 31, 32] that approximates the solution manifold of equations parametrized by control variables $\alpha$ (for example, model parameters) by identifying reduced coordinates for the system, and then constructing equations for these coordinates.

In particular, the reduced basis method (RBM) [33, 34] traditionally works by identifying the reduced coordinates as the amplitudes of a reduced basis informed by previous high fidelity (HF) evaluations of the full system, and creates the reduced equations by projecting the underlying operators using the Petrov-Galerkin scheme [23]. This approach has shown great performance with speed ups of the original calculations usually between 2 and 7 orders of magnitude, yet it has been challenging to successfully apply it when the underlying operators are non-affine [28] or non-linear [25], usually requiring the mitigation of such features by hyperreduction schemes [35] such as the Empirical Interpolation Method [36, 37]. Alternative avenues to construct the reduced equations without explicitly projecting the original full system’ model, but rather by relying on data directly obtained from HF calculations have recently started to be explored [38, 39, 40, 41] (see for example [42] for a related application in other fields). Some of these developments follow a supervised machine learning (SML) approach in where various algorithms are trained to reproduce the dynamics of the reduced coordinates as the control parameters $\alpha$ change. In this work, we explore the application of Genetic Programming (GP) to obtain the reduced equations from data. This is achieved by creating an ensemble of initial possible equations that can reproduce the response of the reduced coordinates to parametric changes, and then performing an evolution algorithm in which the current population of models goes through changes to improve their fitness, measured as their accuracy in capturing the underlying reduced dynamics. We test this framework in two example problems [39] that exhibit non-affine and non-linear characteristics that make them problematic to emulate using the traditional Galerkin projection approach.

II Formalism

II.1 Quantum systems

The first application problem we consider is a modified version of the Gross-Pitaevskii equation [43, 44] (a nonlinear Schrodinger equation used to describe dilute Bose-Einstein condensates):

$$\Big{(}-\frac{d^{2}}{dx^{2}}+\kappa x^{2}+q\rho(x)^{\sigma}\Big{)}\phi^{(i)}(x)=E_{i}\phi^{(i)}(x),$$ (1)

Figure 1: Five wave functions $\phi_{i}(x)$ ($1\leq i\leq 5$) and their associated density self-interaction term $\rho(x)^{\sigma}=\big{(}\sum_{i}^{5}|\phi^{(i)}(x)|^{2}\big{)}^{\sigma}$ for the modified Gross-Pitaevskii equation (1) for $\kappa=1,\ q=2,\ \sigma=2.7$.

in where the Hamiltonian for each wave function $\phi^{(i)}$ consists of a trapping term proportional to $x^{2}$, as well as with a density $\rho(x)=\sum_{i}^{N}|\phi^{(i)}(x)|^{2}$ compromised by all $N$ particles in the system. The control parameters in this case are $\alpha=\{\kappa,q,\sigma\}$. Figure 1 shows the wave functions and their associated density for an example value of the control parameters, while Table 1 shows the ranges of the three parameters explored.

The second application problem we consider is Density Functional Theory (DFT) for a Skyrme type interaction [3]. DFT is used in nuclear physics to describe the atomic nucleus through the mean-field perspective, where the nucleons (protons and neutrons) interact with average densities and currents that are constructed from the nucleon wavefunctions (see  [45] for more details). The Hamiltonian of the system can be written as:

$$\hat{h}_{\alpha}^{(i)}[\Phi]\phi^{(i)}=E_{i}\phi^{(i)},$$ (2)

where the wavefunction of each orbital $\phi_{i}$ interacts with a Hamiltonian $\hat{h}_{\alpha}^{(i)}$ that depends on all of them $\Phi=\{\phi^{(i)}\}_{i=1}^{N}$. This single particle Hamiltonian is derived from the Skyrme effective interaction [46, 47, 48, 49], with the nuclear part is made of time-even densities [50]:

$${\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},$$ (3)

with the subscript $t=(0,1)$ denoting isoscalar and isovector densities, respectively. The parameters of this EDF model the interaction between the wavefunctions and the respective nucleonic densities in question (the kinetic energy density. We usually re-parametrize these constants by terms representing nuclear-matter properties [51, 52]:

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

For this work we focus on ⁴⁸Ca which, within the spherical DFT construction contains 6 independent proton orbitals and 7 independent neutron orbitals. Table 1 shows the ranges of the two varied parameters for this study, $J_{\text{sym}}^{\text{NM}}$ and $L_{\text{sym}}^{\text{NM}}$, which represent the symmetry energy and its slope [52] respectively.

Parameter	Range
$\kappa$	[0.5,3.0]
$q$	[-2.0,2.0]
$\sigma$	[0.5,3.0]
$J_{\text{sym}}^{\text{NM}}$	[20,40]
$L_{\text{sym}}^{\text{NM}}$	[25,65]

Table 1: Parameters varied for each of the two explored problems with their associated ranges. For the DFT case, both the symmetry energy $J_{\text{sym}}^{\text{NM}}$ and its slope $L_{\text{sym}}^{\text{NM}}$ are expressed in MeV. A total of 492 evaluations of the modified Gross-Pitaevskii densities were used for training and 212 for testing. A total of 25 evaluations of the DFT protons and neutrons densities were used for training and 16 for testing.

II.2 Reduced basis method

The first step within the RBM formalism consists of building a reduced basis of the solution to the parametrized equation we are interested in solving:

$$F_{\alpha}[\phi^{(i)}]=0,$$ (5)

where we rewrite the equation in terms of the operator $F_{\alpha}$ which is, for example $F_{\alpha}[\phi^{(i)}]=\big{[}\hat{h}_{\alpha}^{(i)}[\Phi]-E_{i}\big{]}\phi^{(i)}$ in the DFT case 2. For each of the two examples described (1) and (2) we could approximate each wave functions as:

$$\phi^{(i)}(x;\alpha)\approx\hat{\phi}^{(i)}(x;\alpha)=\sum_{k}^{n}a_{k}(\alpha)\phi^{(i)}_{k}(x),$$ (6)

where the few $n$ coefficients $a_{k}(\alpha)$ represent latent reduced coordinates that depend on the control variables $\alpha$, while the reduced basis $\phi_{k}^{(i)}$ is constructed informed by previous HF solutions for different values of $\alpha$. In several recent successful RBM applications to nuclear physics [24, 23, 25, 28, 41, 29], the reduced basis is constructed by performing a principal component analysis [53], or singular value decomposition [54] on the set of HF solutions.

Figure 2: Densities calculated for various values of the controlling parameters $\alpha$ specified in Table 1. The inset plots show the exponential decay of the singular values from the SVD. Panel a) Density $\rho$ for the modified Gross-Pitaevskii Eq. (1). Panel b) Densities of neutrons $\rho_{n}$ and decay of singular values for both protons $\rho_{p}$ and neutrons.

The second step is to obtain the equations that describe the response of the latent variables $a_{k}$ to changes in the parameters $\alpha$. In many RBM applications (including several of the previously cited work), a way to achieve this is to project the equations into a subspace spanned by test functions $\psi_{j}$ [23, 33, 34]:

$$\langle\psi_{j}|F_{\alpha}(\hat{\phi}^{(i)})\rangle=0,\quad\text{for }1\leq j\leq n,$$ (7)

which creates $n$ equations to be solved for the $n$ unknown coefficients $a_{k}$ (in the case of coupled equations, such as the non-linear systems in Eq. (2), there would be $n$ projection equations per original wave-function equation).

In our previous work [23], we successfully constructed a reduced order model for the Gross-Pitaevskii equation (1) in the case of $N=1$ and $\sigma=1$ using the described Galerkin projection scheme (7), yet for bigger powers $\sigma$ for the density, or even fractional powers, this approach proves inadequate to create efficient reduced equations. The challenge lies in the inability to pre-compute the Galerkin projections in the offline stage of building the emulator, since the equations are non-affine and non-linear in $\sigma$ and the wavefunctions. Various terms in realistic EDFs [55, 56, 57] contains fractional powers of densities, which motivates the development of alternative approaches to obtain reduced equations, such as the one we present next.

II.3 Genetic Programming

Figure 3: A plot exemplifying individuals representation in Tree Genetic Programming.

Genetic Programming (GP) is a type of Artificial Intelligence (AI) that falls under the broader category of evolutionary computation (EC), which is a subset of AI inspired by biological evolution mechanisms such as selection, mutation, and crossover to solve problems [58, 59, 60, 61]. Specifically, GP is a population-based stochastic iterative search algorithm in the search space of computer programs. Like other evolutionary meta heuristics, GP evolves a set of candidate solutions (the population) by mimicking the basic principles of Darwinian evolution. The evolutionary process involves an iterative application of a fitness-based selection of the candidate solutions and their variation throughout genetically-inspired operators, such as crossover and mutation [58]. Three of the most typical representations employed in GP are expression trees in tree GP (TGP), linear sequences of instructions in linear GP (LGP) and circuit-type graphs in Cartesian GP (CGP). Here we shall be mostly concerned with TGP, and want to only briefly touch the other two representations. In this work, we will use TGP and will refer to it as simply GP. In GP, the evolving programs are constructed by composing elements belonging to two specific, predefined, sets: a set of primitive functions $F$, which appear as the internal nodes of the trees, a set of terminal input features $T$ and constant values $C$, which represent the leaves of the trees. In the context of SML problem-solving, the trees encode symbolic expressions mapping inputs $X$ to the outputs $y$. Typically, GP is used with the so-called subtree mutation and swap crossover [58]. The latter exchanges two randomly selected subtrees between two different parent individuals. The former randomly selects a subtree in the structure of the parent individual and replaces it with a new, randomly generated tree. Algorithm 1 provides a pseudo-code for the GP, whereas Figure 3 provides a visual representation of GP-models and their mapping between tree and function.

1. Create an $n$-sized initial population of random individuals $P$ 2. Repeat until matching a termination condition, typically the number of generations (a) Calculate fitness $\forall$ individual $i$ in $P$ (b) Create an empty population $P^{\prime}$ (c) Repeat until $P^{\prime}$ contains $n$ individuals i. With probability $P_{c}$ choose the crossover operator. Alternatively, mutation will be selected ($1-P_{c}$) ii. Select 2 individuals using some selection algorithm iii. Apply the operator selected in 2.c.i to the individuals selected in 2.c.ii iv. Insert the individuals from 2.c.iii into $P^{\prime}$ (d) Replace $P$ with $P^{\prime}$ (and copy the elite(s) from $P$, if elitism is used) 3. Return the best individual in $P$

Algorithm 1 Pseudo-code for the typical GP algorithm.

In an ML context, GP provides a user with a completely different experience in terms of transparency and interpretability of the final models when compared to more “black-box” models of, say, deep artificial neural networks (ANNs). This is because GP typically outputs symbolic expressions or tree structures representing programs, which are generally interpretable and can be analyzed by humans. By manipulating discrete structures, the inner workings of a GP-based model are clearer, and the relationships between features and predictions are explicitly defined. As a result, evolved solutions can offer valuable insights into how the model arrives at its predictions, allowing humans to understand the decision-making process and relate this to their domain-knowledge.

Although GP is not the only “white-box” approach in the spectrum of ML methods, it allows for unprecedented flexibility of representations and modularity, making it suitable for numerous tasks: classification [62, 63, 64], regression [65, 66, 67], feature engineering [68, 69], manifold learning [70], active learning [71], image classification [72, 73, 74, 75], image segmentation [76], image enhancement [77, 78], automatic generation of ML pipelines [79] and even neural architecture search [80, 81]. The area of explainable AI receives consistently more attention from both practitioners and researchers [82, 83] and GP has gained popularity where human-interpretable solutions are paramount. Real-world examples include medical image segmentation [84, 85], prediction of human oral bioavailability of drugs [86], skin cancer classification from lesion images [87], and even conception of models of visual perception [88], etc. Besides being able to provide interpretable models, there is evidence that GP can also help to unlock the behaviour of black-box models [89, 76].

We proceed to apply GP to learn equations for the reduced coordinates the particle density:

$$\rho(x;\alpha)=\sum_{i}^{N}|\phi^{(i)}(x)|^{2}\approx\hat{\rho}(x;\alpha)=\sum_{k}^{n}a_{k}(\alpha)\rho_{k}(x).$$ (8)

That is, for both problems in Eq. (1) and Eq. (2), we learn equations that can describe how each coefficient $a_{k}$ changes as we change their respective parameters $\alpha$ within the defined ranges in Table 1.

In both cases we select a total of three basis $n=3$. For the DFT case we analyze the proton density $\rho_{p}$ as well as the neutron density $\rho_{n}$. The reduced basis $\rho_{k}(x)$ for each of the respective densities is obtained following a singular value decomposition on a set of HF solutions for the training parameters (see Table 1). Note that, since the original equations Eq. (1) and Eq. (2) can’t directly be solved by only considering the density (they involve all wave functions), this GP approach is effectively learning reduced equations that would not be possible to obtain through a direct Galerkin projection.

Table 2 lists the hyper-parameters (HPs) used for GP, along with cross-validation settings. The HPs were selected following common practices found across the literature to avoid a computationally demanding tuning phase. $R^{2}$ was used as a fitness function [90, 91, 92] as it was found to converge faster, generalize better, even when only a few data points are available. A common practice in GP is to protect operators to avoid undefined mathematical behaviour of the function set by defining some ad-hoc behaviour at those points like, for instance, returning the value 1 in the case of a division by zero to make it possible for genetic programming to synthesize constants by using $x/x$ [58]. However, it was shown that these techniques present several shortcomings in the vicinity of mathematical singularities [90]. In this study, programs that produced invalid values were automatically assigned a bad fitness value, making them unlikely to be selected. By diminishing the selection attractiveness of solutions that produced invalid values, we expect to obtain models whose fitness landscape is less sharp [93].

For GP we used the commercially available DataModeler system [94]. The parameter set used is shown in Table 2.

Parameters	Values
№runs	16
Population size	300
Functions ($F$)	{+, -, x, /  -x, CF, $x^{2}$, $x^{-1}$, RP, $\sqrt{x}$, $e^{x}$, $log(x)$, $x^{n}$, $x^{1/3}$, $x^{3}$}
Selection	Pareto tournament selection size 30
Genetic operators	{subtree crossover, subtree mutation, depth preserving subtree mutation}
$P(C)$	0.9
$P(M)$	0.1
Maximum complexity	1000 (Visitation Length [95])
Stopping criteria	5 minutes

Table 2: Summary of the GP hyper-parameters. Note that $P(C)$ and $P(M)$ indicate the crossover and the mutation probabilities, respectively. CF represents the continued fraction operator $d/(b+c/a)$ and RP represent the rational polynomial function $(b+d+f)/(a+c+e)$.

To contrast the GP ability to learn the equations from data, we also implemented a linear regression model. In the context of SML, linear models are a class of algorithms that assume a linear relationship between the input variables (features) and the output variable (target), in this case the controlling parameters $\alpha$ and the reduced coefficients $a_{k}$, respectively. To make a prediction, the input features are multiplied by the learnable weights, summed together and added a learnable constant called the intercept or bias [96]. Formally: $y=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\ldots+\beta_{n}x_{n}$, where $y$ is the target, $\beta_{0}$ is the intercept, and $\beta_{k}$ is a learnable weight associated with input feature $x_{k}$. The learnable weights and the intercept are usually estimated from the data using Ordinary Least Squares (OLS) estimation procedure to minimize the sum of squared residuals (i.e., difference between the observed target values in the training dataset, and the values predicted by the linear model). Formally: $\min_{\beta}\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}$, where $\hat{y}_{i}$ is the model’s prediction for data instance $i$.

Figure 4: Performance analysis of the SML approaches for the modified Gross-Pitaevskii Eq. (1). Panel a): Computational Accuracy vs Time plot [28] that shows for each of the 212 testing evaluations the norm of the residual in the vertical axis and the calculation time in the horizontal axis for cubic linear regression with L1 regularization (dark red), with L2 regularization (dark green), and with GP (blue). The corresponding dashed vertical lines are the median time for each method, while the dotted horizontal lines the median error. The horizontal dashed black line shows the median error obtained when optimizing the coefficients in Eq. (8) directly against the test solution. The vertical grey solid line shows the median calculation time of the HF solver Python implementation, which is more than 4 orders of magnitude slower than any of the three methods. Panels b) and c) show the evolution of the first coefficient of the reduced basis $a_{1}$ as the parameters $\kappa$ and $\sigma$ change, respectively, while the other parameters are held constant. The black dashed line shows the optimal value of the coefficient as the parameters change. The dashed vertical grey line shows the boundary on which the SML methods were trained, after which the calculations are effectively extrapolating.

In order to constraint the complexity of the model and improve generalization ability, regularization terms can be added to the sum of squared residuals. L1 regularization promotes models’ sparsity and is useful to tackle high dimensional tasks (i.e. with many input features), especially if many are expected to be irrelevant. L2 regularization is frequently used when the input features are expected to be informative but their impact has to be reduced to avoid overfitting, and/or are highly correlated (multicollinear). Formally, the L1 term is given by $\lambda_{1}\sum_{k=0}^{n}|\beta_{k}|$, whereas L2 term is given by $\lambda_{2}\sum_{k=0}^{n}\beta_{k}^{2}$. For creating high-order interactions between the original input features $X$, one can generate a new feature set $X^{\prime}$ consisting of all polynomial combinations of the features up to a given degree. Provided $X^{\prime}$, a linear model can learn non-linear relationships. In this study, we use regularized linear regression with polynomial terms and we rely on scikit-learn implementation [97, 98].

Note that regularized linear regression has been long used in scientific applications. Sparse-identification of non-linear dynamics (SINDy) is a popular method for discovering governing equations in dynamic systems [31]. Essentially it utilizes linear regression with the L1 regularizing term to promote sparsity, yielding smaller models, which can offer better interpretability, which is considered an appealing property. SINDy has been successfully applied to find governing equations in many different applications [99, 100, 101, 102, 103]. Nevertheless, SINDy-like approaches present a conspicuous disadvantage: all of the non-linear relationships to be included in the model search must be manually pre-defined, which requires some degree of expertise. This can be particularly limiting when applied to a problem where valuable non-linear relationships are unexpected, and is one of the advantages we expect from GP.

III Results

Figure 5: Same as Panel a) of Figure 4 but for the case of the neutron density $\rho_{n}(x)$ in DFT for 16 values of $\alpha$ in the testing set. The vertical grey solid line represent the typical spherical DFT calculation time in fortran. This time is faster than for the simpler case in Figure 4 since research codes, like the one used here for DFT, have been optimized for performance. Yet, all SML models which show similar times to those obtained in the case of the modified Gross-Pitaevskii equation are faster by almost 4 orders of magnitude. The horizontal dashed black line shows the median of the optimal coefficients. We observed very similar results for the proton density $\rho_{p}(x)$.

Once the models were trained, we tested their generalization ability on the previously unseen combinations of parameters belonging to the respective test sets. Given a set of parameters $\alpha_{i}$ from the test set, we quantify the accuracy performance as the root mean squared error (or $L^{2}$ norm) between the density obtained from a high-fidelity solver (HF) and the one obtained from the reduced coordinates equation $\phi(x;\alpha_{i})\approx\sum_{k}^{n}\hat{a_{k}}(\alpha_{i})\phi_{k}(x)$, where $\hat{a_{k}}$ were predicted by the respective SML models as a function of $\alpha_{i}$. We quantify the speed performance as the time it takes for each procedure to estimate the reduced coordinates given the target $\alpha_{i}$. Figure 4 shows the performance of each SML model (linear regression and GP) for the modified Gross-Pitaevskii equation. Panel a) shows the accuracy and speed of the three methods. The horizontal dashed black line shows the median error obtained when the coefficients in Eq. (8) are optimized to directly reproduce the test solutions, which serves as a lower bound for the error, only possible to be improved by increasing the number of kept principal components beyond $n=3$. Both GP and LR3-L2 attain residuals comparable with this bound. In terms of speed, all three methods are more than 4 orders of magnitude faster than the HF solver.

Panels b) and c) of Figure 4 show how the first coefficient of the reduced basis $a_{1}$ evolves as the parameters $\kappa$ and $\sigma$ change, respectively, while the other parameters are held constant. This is known as ceteris paribus or isolated effect and is frequently used when interpreting the coefficients of a linear regression model, as it assesses the effect of a particular input feature (predictor) on the response (dependent variable) while keeping all other predictors fixed. Analogous to panel a), the black dashed line represents the evolution of the optimal value of $a_{1}$ as the parameters $\alpha$ change. The GP has appreciably better performance in reproducing this optimal coefficient that the other two methods, and it excells particularly in extrapolating beyond the dashed vertical line which denotes the maximum value of the parameters to which all the models were trained. This suggests that the GP has learned the reduced dynamics with less overfitting.

Figure 5 shows analogous results as panel a) of Figure 4 for the DFT neutron density $\rho_{n}(x)$ of ⁴⁸Ca (comparable results were obtained for the proton density). Similar to the performance shown in the case of the modified Gross-Pitaevskii equation, the GP method (as well as LR3-L2) is able to learn the dynamics almost perfectly, showing accuracy that is almost on top of the optimized coefficients indicated by the black dashed line. The computation times of all SML methods are comparable to the easier problem in Figure 4. Indeed, it is remarkable that GP is able to find effective equations that explain the change in the nuclear densities as the two parameters $J_{\text{sym}}^{\text{NM}}$ and $L_{\text{sym}}^{\text{NM}}$ are widely varied. The reduced order model we have constructed not only is agnostic to the high dimensional grid-space in the coordinate $x$, but it is also learning the solution manifold of the densities without having to directly track the underlying wave functions, and in the case of DFT the other densities involved. These results poise GP as a promising tool for drastically reducing the computation cost of existing complex nuclear calculations, as well as for model discovery directly from experimental data.