Trees and Forests in Nuclear Physics

2020·Arxiv

Abstract

Abstract

We present a simple introduction to the decision tree algorithm using some examples from nuclear physics. We show how to improve the accuracy of the classical liquid drop nuclear mass model by performing Feature Engineering with a decision tree. Finally, we apply the method to the Duflo-Zuker model showing that, despite their simplicity, decision trees are capable of improving the description of nuclear masses using a limited number of free parameters.

I. INTRODUCTION

In recent years, there has been an enormous growth of new statistical tools for data science [1, 2]. Although these methods are extremely powerful to understand complex data and detect novel patterns, they are still rarely adopted by the nuclear physics community. Only a few groups are currently pioneering the applications of these methods to the field. These topics have been recently discussed in a series of workshops on Information and Statistics in Nuclear Experiment and Theory (ISNET). Recent developments in this field are documented in the associated focus issue published in Journal of Physics G [3]. The aim of this guide is to illustrate an algorithm used widely in data analysis. Similarly to our previous guide on bootstrap techniques [4], we present the decision tree starting from very basic models, then finally apply it to more realistic problems, like improving models for nuclear mass predictions.

Decision trees are already implemented within major experimental collaborations, such as MiniBooNE, to improve the performances of particle detectors [5, 6], but they are not yet widely used in low energy nuclear physics, where they could help to analyse both experimental data [7] and theoretical models.

Following the notation and terminology of Leo Breiman’s paper Statistical Modeling: The Two Cultures [8], we want to investigate a process f that transforms some input X into an output Y . That is to say, f is a function:

where the input X can be quite general, from images to a table of data, while the output Y can be a discrete or continuous set. In the first case we speak of a classification problem, in the latter of a regression problem.

Rather than focusing on investigating the fine details of the process f with many restrictive assumptions (an approach that is named data model culture in [8]), we consider f as a black box mapping X to Y , and we try to approximate it. That is, we give up trying to investigate all the fine details of f and we focus on finding a representation (or approximation) ˜f for f. ˜f is called a model and it is a function with the same domain X of the process f and codomain ˆY :

The process ˜f depends on variables (usually named features), parameters (coefficients that can be learned with the algorithm) and hyper-parameters, that are set before training the model (and thus are not learned). We will present an extended discussion on how to select the features of the model to improve performances in Sec. II C. Another goal for the feature selection process is to reach a representation as parsimonious as possible.

Since “all models are wrong, but some are useful” [9], it is necessary to introduce a definition of what a good model looks like in order to pick the best one out of a set of possible candidates. Or in other terms, we need to assess how faithfully ˜f represents the process f. Mathematically, the goal for training a model ˜f is to minimise a particular scoring function, sometimes improperly called “a metric”. Without loss of generality, we are considering only the minimisation problem: changing the sign of a scoring function to be maximised reduces the problem to a minimisation task.

For example, a natural choice for the scoring function is the mean squared error (MSE) or variance, that is the difference between the predicted value ( ˆY ) and the observed, experimental data (Y ) [10]:

It is worth noting that this is not the only option, and within the Machine Learning literature we encounter other scoring functions as the logarithmic mean squared error (MSLE):

or the median absolute error [10]:

Different scoring functions correspond to different modelling choices and the importance we assign to specific subsets of the database. The use of MedAE would be more appropriate to obtain a model that is robust to outliers: a few poorly described experimental points will not alter significantly the performances. In the current work, we have chosen the mean standard error (or equivalently its square root, RMS) which is the default in most libraries. Given the high accuracy of measurements in nuclear physics, especially for masses as discussed here, we do not need to worry about possible outliers in our data sets and MSE therefore represents a reasonable choice.

Another important aspect in building a model is the decision on the tradeoff required between model performances and explainability. That is, the choice between (possibly) better performances with the chosen scoring and easier explanations of the model in plain language. Among the regressors usually considered to be explainable are linear regression and decision trees. However, some recent research allows explaining (although approximately) even the results from algorithms deemed black-box, such as neural networks, or such as gradient boosting in explainable models like linear regression [11] and simple decision trees [12].

In this guide, we chose to illustrate decision trees because they retain explainability, they do not rely on the assumption of linearity nor on the linear independence of the features, and they are not significantly affected by monotonic transformations (no input data scaling is required, nor monotonic transformations like taking the logarithm or the square of one variable). Also, decision trees are the key elements in building other regressors like Random Forests [13] or Xgboost [14] that usually perform better with regard to scoring.

Last but not least, an important aspect of modelling is the balance between the complexity of the chosen model and the generality of the results. As an analogy, it is useful to consider the problem of approximating N experimental distinct points using a polynomial. A complex polynomial of degree N will be able to describe perfectly the data. Whenever some new data are added, the perfect description will (in general) no longer be true. A correct assessment of the performance of a regressor should be performed on unseen data, i.e. data that were not used during the training.

A common practice to estimate the performance on unseen data is the k-fold cross validation, with . In essence, the data are permuted and then separated in sets of size k with each subset (fold) roughly of the same cardinality. The model is then trained on 1 subsets and validated on the subset not used while training. As an extreme case, when 1, all data but one are used for training and the performances are assessed on only one datum. This scheme is called “leave-one out validation” [15].

In the following sections, we will illustrate the behaviour of decision trees using some nuclear mass models. The article is organised as follows: in Sec.II, we provide an introduction to what a decision tree is, using very simple examples. In Sec.III, we introduce the nuclear models to which we will apply the decision trees. The results of our work are presented in Sec.V and we illustrate our conclusions in Sec.V. In the Supplementary Material, we provide the Python script used to perform the calculations. The script has been structured in the same way as the current guide to facilitate its usage [43].

II. DECISION TREE

With a decision tree, the function is approximated with a step function with n steps as

with Ωwhere d is the number of features, ) is the indicator function:

and Ωare half-planes in R.

Any measurable function can be approximated in terms of step functions [16], thus the approximation is justified as long as the function f is expected to be measurable. That is, using enough step functions we can approximate any measurable function.

Each step function required to build the model ˜f (the tree) is called a leaf, thus the number of leaves of the model corresponds to the number of step functions employed. In order to determine the optimal values for the and Ωof Eq.6, one should provide a splitting criterion; for example, being an extreme value (maximum or minimum) for a given function L. Here, we decide to minimise the norm of the difference between f and ˜f, that is:

This function should be chosen to approximate the scoring function selected; then determining the extremes of L guarantees that we have optimised the desired scoring function. Since in this guide we chose as a scoring function the MSE, a natural choice for the splitting criterion is the norm of Eq.8; we will use it through all the following examples.

We are going to focus on the CART algorithm as presented in [1, 17]. Calculating all the possible splits for all the features to get the optimal ˜f as in Eq.6 is computationally unfeasible for large data sets. For this reason, a greedy approximation is used for training the decision tree: at every iteration of the algorithm, a local optimal split is selected. This is a heuristic approach and there is no guarantee of converging to the global optimum.

At the first step of the algorithm, all the possible splits for all the features are explored, and the best split (that minimises L) is selected. Then, all the data are split between the two leaves (a leaf for each half-plane). Then, for every leaf, the procedure is iterated until a stopping criterion is reached. There are many different stopping criteria: a leaf can not be further split if it contains only one observation or if all the features are constant. However, the training is usually stopped as a modelling choice to avoid poor performance on unseen data (overfitting): once a given number of leaves or a maximum depth (that is, a maximum number of splits between the root and any leaf) are reached, the algorithm stops. Alternatively, leaves are not split if they contain fewer than a specified number of observations. In this process, some of the features may have never been used; in this case, they are irrelevant to the model, and their absence from the input will make no difference.

In the next subsections, we provide some examples of how a decision tree operates, by showing artificially simple examples of trees with few features and very few leaves that can be easily understood. The more realistic cases will be illustrated in Sec.III.

A. A single variable example

As a first example, we illustrate the first iteration of a decision tree, splitting the data (a single feature) into only two leaves. We build an example training set defined as the union of two sets and contains random points uniformly distributed in [0, 0.5] and analogously with points in (0.5, 1]. Notice that . In this case, there is only one feature, so d = 1. In this example as well as in the following one, we will directly code the necessary steps to obtain the desired solution; the more advanced reader can skip these two cases to Sec.II C where an existing library is used.

For the images through f of ), here named ), we use a Gaussian distribution with mean 0 (1) and a variance of 0.05 (0.1). The training set is illustrated in Fig.1. All figures presented in this article have been realised with the Python [18] library matplotlib [19].

By construction, the data set can be fully described using a decision tree with only two leaves, that is:

By visually inspecting the data shown in Fig.1, we notice that the current data belong to two groups and a single split along the x axis will be enough to describe them. In more advanced examples, the number of leaves will be selected algorithmically. To train a decision tree means to determine the function given in Eq.9, in such a way that

FIG. 1: The full dots represent the original data set, while the line represents the approximating function ˜f. See text for details.

FIG. 2: Evolution of the L norm defined in Eq.10 as a function of the splitting point . See text for details.

is minimal. L is equivalent to Eq.3 apart from a global scaling factor. The sets Ωand Ωare defined as

It is easy to prove that the constant which best approximates (in terms of norm) a set of values is the average of

the values, ¯x.

The optimal value for (0.5) is obvious when the generating process for the data is known, but how do we determine it when f is not known? The answer is reported in Fig.2, where we plot L as a function of . The optimal value is the one that minimises L.

For this particular case we obtained = 0.493, = 0.003 and = 1.007. More details for reproducing the results are provided in the Supplementary Material. Thus the decision tree reads:

FIG. 3: Graphical representation of the data set X defined by Eq.12. The squares correspond to Y = 10, the dots Y = 1 and the triangles to Y = 0. The solid line corresponds to the first optimal split, the dotted line to the second split performed by the decision tree. See text for details.

Following these simple steps, we have built a model ˜f that is able to provide a reasonable description of the main structure of the data, i.e. we recognise that the data are separated in two groups. This is represented by the solid line in Fig.1. We say that this tree has two leaves since we have separated the data into two subgroups. Notice that the MSE is not exactly equal to zero since there was some noise in the generated data.

B. A two variables example

We now consider a slightly more complex problem with two variables . The aim of this example is to familiarise with the concept of multiple splits to treat a complex problem via simple operations. As in the previous case, we will apply the basic steps to explicitly build a decision tree. We generate a new set of data points as:

with response:

Graphically, we represent this data set in Fig.3. The data are clearly clustered (by construction) in three regions of the plane. The aim of the current example is to illustrate how to perform successive splits to correctly identify these regions.

We apply a two-step procedure: firstly, we separate the data along the direction. Following the procedure highlighted in the previous example, we create a model in the direction of the form

We now perform a systematic calculation of the L norm looking for the value that leads to its minimal value. We refer to the Supplementary Material for details. We find = 0.482 as the value for dividing the plane. By observing the data, we see that there is no gain by adding further splits in this direction. We will come back to this aspect in the following sections. We can further refine the model by adding an additional separation along the direction. The procedure follows the same steps as before and we find that = 0.495. The result is reported in Fig.3.

An important quantity for analysing the model and for assessing the importance of its input variables is an estimate of the feature importance. This is particularly relevant for the ensemble methods that rely on decision trees: while with a single decision tree the role of the features is obvious once the tree is visually represented (see for example Fig. 4), it is unpractical to represent all of the decision trees in a Random Forest (there may be thousands). Also, a feature may participate multiple times in different trees, so a definition of importance should take this into account.

Starting from Fig. 4, we want to assess the importance of the features to the building of the decision tree. We recall here that by features we mean the variable of the model. In this particular example, we have considered as a natural choice, but one could also consider other combinations: . We refer the reader to Sec.II D for a more detailed discussion.

Following Ref. [10], we calculate the feature importance in the following way: for each split s in the tree, we calculate the weighted impurity decrease as

N is the total number of observations, is the number of observations at the current node, is the number of samples in the left leaf, and is the number of samples in the right leaf. I represents the impurity (in our case, MSE), with the subscripts having the same meaning as before.

By inspecting Fig.4, we observe that for the first split there are in the current node (the root of the tree) as many observations as the total, that is = 100. The initial impurity (MSE) is 22.782, = 50, right impurity is 0, left impurity 0.25. Thus we obtain

For the second split, we get:

Normalising the total weighted variation to 1, we obtain that the column has importance equal to 99.5% for the model, while has an importance of 0.5%. If the variables entered in different splits, the relative importance would be summed.

Estimating feature importance is fundamental for improving the quality of the model: by discarding irrelevant features, i.e. features that are not reducing the impurity in the tree, more parsimonious models can be trained. This is especially useful for models involving hundreds (or thousands) of features. For example, anticipating the results of the following section, we see that in Figure 10, a simpler model could be obtained using only 6 features instead of 9. Whenever features are generated for improving the model, a critical assessment of their relevance should be performed.

After these examples that were easily implemented with a few lines of code, in the next sections (and for realistic problems, in terms of the number of features and the number of possible leaves) we are going to rely on existing libraries.

C. A two variables example (revisited)

In this section, we make use of the function DecisionTreeRegressor from the Python package scikit-learn to determine the structure of the tree, using all default hyper-parameters apart from the number of leaves. For the sake of simplicity, we still consider as features of the model and we also impose the number of leaves to be three. In more advanced examples, we will let it be a free parameter.

In this case, contrary to the previous example, we do not need to decide if the split along the direction happens before the one along the or vice-versa: all possible splits on the available data for all features are explored with the algorithm.

The algorithm used to perform such a split can be represented as in Fig.4 using the scikit-learn package [10], but an analogous result could have been obtained using R [20], and the libraries rpart [21] for model training and rpart.plot [20] for visualisation.

The visualization of a scikit-learn tree consist of a series of boxes counting basic information: the value of the variable at which the separation takes place, but also the mean value of the data (named value) and the impurity (MSE in our example). It also provides information concerning the amount of data grouped in each leaf. In this case, the tree has a total of three leaves. As the MSE is zero (thus, minimal) on each leaf, adding extra splits to the model would not lead to any real gain in the description of the data, but it would only increase the model complexity.

FIG. 4: Decision tree with two variables obtained using the scikit-learn package [10] for the data set reported in Fig.3.

TABLE I: Evolution of MSE with the number of leaves of the tree for the data set shown in Fig.5.

D. The importance of Feature Engineering

In the previous examples, we have approximated the data using simple step functions; although this choice is mathematically justified, the problem is that the approximation may lead to a single observation per leaf, with the result that the generalisation on unseen data may be unsatisfactory.

To overcome the problem, and possibly to make the models easier to explain, it is important to explore the data and apply convenient transformations on the input variable for the model that may highlight some patterns. We consider as an example the case of a two variable data-set obtained as follow:

where the response is chosen as

In the left panel of Fig.5, we illustrate the data points using Cartesian coordinates. By using DecisionTreeRegressor, we build a series of decision trees as a function of the number of leaves. In Tab.I, we report the MSE of the tree for various number of leaves.

From the table, we see that the lowest MSE is given by a 7-leaves tree. This tree is quite complex, with leaves containing only 2 or 4 observations. For illustrative purposes let’s consider the case with 3 leaves, which corresponds to MSE = 0.073. The splits are performed along the direction at = 0.742 and along the direction at = 0.979. The result is reported using solid and dashed lines in the left panel of Fig.5.

By inspecting the data, we realise that by applying a unitary transformation from Cartesian to polar coordinates:

we can highlight a very specific structure in the data. The result of such a transformation is illustrated in the right panel of Fig.5. Using these new variables in the model obviously helps the performances: using a single split, i.e. a

FIG. 5: Graphical representation of the data set (X,Y): on the left panel the Cartesian coordinates are used while on the right panel the data are represented in polar form. The triangles correspond to Y = 1 and the dots to Y = 0. The solid lines represent the cuts done to reproduce the data.

tree with two leaves, we obtain a total MSE of zero using fewer leaves and with features that are easier to understand. Only one feature is relevant (), so while the possible splits on the other features are explored with the algorithm, they are irrelevant and do not play any role in the model. Notice that by construction, there is only radial dependence and no dependency on the phase. Thus we obtain a more parsimonious model by applying Features Engineering.

E. Random forests

We conclude this section on decision trees by introducing the concept of a Random Forest (RF). Following [1, 13, 22], we define a RF as an ensemble of decision trees. The first step for training a RF is to use bagging, also named bootstrap aggregating [23, 24]. Given a training set , a uniform sampling with replacement is performed, obtaining two data sets: one containing the sampled data (with repetitions), , and one containing the data that were never sampled, named out-of-the-box (OBB), contains roughly one-third of the initial observations.

In fact, given a number N of observation, assuming no repetitions in the data, the probability for a datum of not being extracted at the i-th draw is simply = 1 . Thus the probability of not being extracted after N draws, having replacement and assuming all the draws to be independent, is:

The second step of the RF algorithms consists of selecting a random subset of the features of . The number of features used is an adjustable hyper-parameter of the algorithm; with [10] for example, the user provides the fraction in (0, 1) out of the total number of features that will be used for each tree. A decision tree is then built on using only the selected features, and the performances are estimated on . In [10] the procedure is fully automatic and many trees can be trained in parallel, but the estimation on is not performed by default (but the option can be activated). Repeating the bagging and the tree training for a number T of trees (another adjustable hyper-parameter) and averaging the response ˆY over the ensemble of predictions, a Random Forest is obtained.

The bagging and the random selection of features are tools to inject noise in the training data. The noise can be reduced by averaging on the response of each tree, and this empirically improves the performance of the regressor [13, 23]. For the theoretical reasons why the performances are better after injecting some noise in the system, the reference is [13].

Intuitively, the trees in the forest should not all provide the same response, otherwise averaging on all the trees would be of no benefit. Consider for example the dataset X of Eq. 12: if all the trees are trained on the same data with the same input features, then they will all provide the same output, and the random forest will be equivalent to a single decision tree.

A reason that made Random Forests very popular is that they can be trained on data sets with more features than observations without prior feature selection, a characteristic that made then a relevant tool, for example, for gene expression problems in Bioinformatics. For more details, see [25].

III. NUCLEAR MASS MODELS

Before applying the decision tree to a more realistic case, we now introduce the nuclear models we are going to study. According to the most recent nuclear mass table [26], more than 2000 nuclei have been observed experimentally. To interpret such a quantity of data, several mass models have been developed over the years [27] with various levels of accuracy. For this guide, we selected two simple ones: the Bethe-Weizs¨acker mass formula [28] based on the liquid drop (LD) approximation and the Duflo-Zuker [29] mass model. The reason of this choice is twofold: the models contain quite different physical knowledge about the data, for example, the lack of shell effects in LD case, but they are relatively simple and not CPU intensive, thus giving us the opportunity to focus more on the statistical aspects of the current guide.

A. Liquid drop

Within the LD model, the binding energy (B) of a nucleus is calculated as a sum of five terms as a function of the total number of neutrons (N) and protons (Z) as

where A = N + Z. The set of optimal parameters have been tuned in Ref. [4]. These parameters are named volume (), surface (), Coulomb (), asymmetry () and pairing () and they refer to specific physical properties of the underlying nuclear system [28].

B. Duflo-Zuker

The Duflo-Zuker [29] is a macroscopic mass model based on a generalised LD plus the shell-model monopole Hamiltonian and it is used to obtain the binding energies of nuclei along the whole nuclear chart with quite a remarkable accuracy. The nuclear binding energy for a given nucleus is written as a sum of ten terms as

We defined 2and . The ten different contributions can be grouped into two

categories: in the first one we find terms similar to the LD model as Coulomb (), symmetry energy () and pairing . The other parameters originate from the averaging of shell-model Hamiltonian. See [30] for more details. The model described in Eq.21 is usually referred as DZ10 and its parameters have been recently tuned in [31]. Within the literature, it is also possible to find other versions with extra parameters [32], but we will not consider them here for the sake of simplicity.

In Fig.6, we illustrate the behaviour of the residuals E(N, Z) obtained with the two mass models i.e. the difference between the empirical data and the models predictions. We assume that nuclear data [26] have negligible experimental error compared to the model, and we discard all data having an uncertainty larger than 100 keV. This is a reasonable assumption to be made since the typical discrepancy between models and data is usually one or two orders of magnitude larger than the experimental errors [27]. See discussion in [4] for more details.

In each panel of Fig.6, we also provide the root mean square (RMS) deviation . We thus see that DZ10 is roughly one order of magnitude more accurate in reproducing data than the simple LD. The detailed analysis of these residuals has been already performed in [4, 31] showing that they do not have the form of a simple white noise, but they contain a degree of correlation.

FIG. 6: Residuals (expressed in MeV) obtained using the liquid drop model (left panel) and the DZ10 model (right panel).

IV. RESULTS

We apply the decision tree to the case of nuclear data. Using the same notation adopted in the previous examples, the input X is now a matrix with 3 columns N, Z, A while the response ˆY is the residual.

As specified before, the goal is to minimise the RMS on unseen data or, in other words, learning without overfitting. While it appears obvious that a tree with only one leaf, which means replacing all the values of Y with the average Y , or with as many leaves as there are observations are not very useful, determining the optimal value for the number of leaves is not straightforward. The approach is empirical: experimenting with a reasonable set of values for the number of leaves, and pick the best results according to the cross-validation.

With only one adjustable hyper-parameter like the maximum number of leaves, exploring the parameter space is straightforward: all the possible values are tested, with a cost of M cross-validated models, where M is the number of possible values for the number of leaves. In our example, exploring trees with a maximum number of leaves between 2 and 500 implies cross-validating for M = 499 models.

On the other hand, with regressors with many adjustable parameters, as for example Xgboost [14], exploring the hyper-parameter space is more challenging. For example, with 10 hyper-parameters, exploring M values for each of them means exploring a grid with 10points. In this case, it is better to use dedicated libraries [33].

As a first application, we train a simple decision tree over the residuals of the LD model as shown in the left panel of Fig.6. In Fig.7, we illustrate the evolution of the MSE as a function of the number of leaves. For sake of clarity we truncated the figure to 50 leaves, the full plot can be found in the Supplementary Material.

From Fig.7, we notice that the optimal number of leaves is four. The structure of the tree is reported in Fig.8. By inspecting the splits of the data, we notice that the main feature of the data is associated with the neutron number N. The tree splits the nuclei in super-heavy (A > 219) and non-super-heavy. Then it further splits into very neutron-rich and not. Finally, the tree separates out the remaining nuclei into two groups: light and heavy. Having the optimal tree, we now translate it into a simple code. Here we use Fortran, but any other computer language can be used with no difficulty.

FIG. 7: Evolution of the MSE as a function of the number of leaves. The dot corresponds to the absolute minimum. See text for details.

FIG. 8: Decision tree for LD model using only three features N, Z, A. See text for details.

Using the previous code, we now calculate the nuclear binding energies as

where represents the binding energy calculated with the decision tree. By comparing the predicted masses obtained with Eq.22 with the experimental ones, we obtain an RMS of = 2.467 MeV. This is what is called training error, which is the RMS between the response and the predictions of the model trained on all data. A more conservative estimation that should be preferred is the validation error on unseen data, i.e the RMS estimated on data that were not used during the training. In this case, = 2.925 MeV.

FIG. 9: Improved decision tree for LD model using Feature Engineering. For the sake of clarity and readability, the impurity (MSE) was omitted.

It is possible to further improve on this result, by using Feature Engineering as discussed previously. To this respect, we provide some additional information to the tree: . . . . The full list of features can be seen from Fig.10. By inspecting Eq.20, we observe that these features are already used to build the LD model and as such we help the decision tree to identify new patterns in the data. It is worth noting here that other features may be used instead, but a monotonic transformation of existing features (like if we are using A) will provide little to no performance improvement. See for example [34] for an empirical discussion of the topic. Identifying patterns into the data is of great help since it may lead (in complex cases) to better solutions.

In Fig.9, we report the structure of this new tree. The optimal number of leaves is 9. By implementing this tree into a simple numerical code, as done previously, and applying it to the LD residuals we obtain a slight improvement. The total RMS over the entire nuclear chart now falls to = 2.069 MeV (on unseen data, = 2.881 MeV).

Although the decision tree performs less well (in terms of RMS) than a more complex neural network [35], we can still use it to identify possible trends in the data set. By inspecting Fig.10, we observe that not all the 9 features have been used to build the code. In Fig.10, we illustrate the relative importance of the features of the LD model, calculated using Eq.14. We see that the proton fraction Z/A is more important than the individual number of neutrons and protons. It is interesting to note that the decision tree is not affected by even/odd nuclei. This may imply that either the simple pairing term in Eq.20 is enough to grasp the odd-even staggering, or the granularity required to the tree is too high, leading to a number of leaves comparable with the number of data points, or other features can surrogate the odd-even features. We also observe that in this tree the total number of nucleons A and the proton fraction Z/A are as important or more than the number of neutrons N. This clearly explains why the performances of this new tree have improved compared to the one given in Fig.8. A detailed understanding of the trend in the data would require a more in-depth analysis, and so we leave it for future investigations.

In Fig.11, we present a graphical illustration of the energy corrections found by the decision tree for the different nuclei along the Segr´e chart. This figure is an alternative way to represent the various leaves shown in Fig.9. We observe that we have 6 major splits along the valley of stability where we find light, medium-heavy and heavy nuclei. The latter are then still separated into 4 smaller groups. The other cuts occur along the region of proton-rich and neutron-rich, thus the edges of the chart. From this general overview, we may conclude that the residuals of the LD model are quite homogeneous (only two separations) along the valley of stability up to medium-heavy nuclei. Outside

FIG. 10: Relative importance of the features (reduction in impurity normalised) in the Liquid Drop model. The features N/A, (equal to 1 if N or Z is even and 0 otherwise) were not used in the model and as a consequence they have zero importance.

FIG. 11: (Colours online) Graphical representation of the splits done by the decision tree illustrated in Fig.9 on the Segr´e chart of nuclei. The various zones correspond to the energy corrections expressed in MeV derived from the decision tree to the LD model as a function of N, Z.

this range, the number of splits increase since the tree identifies a larger variation in the data. This may imply some missing physics in the model (choice of features) for these particular regions of the chart.

Having seen how the decision tree works for a schematic model as LD, we now apply it to the more sophisticated DZ10. We adopt the same features as shown in Fig.10 to obtain the best performances. Since the structure of the residuals is different there is no a priori reason to use such features, but for the sake of simplicity of the current guide, we keep them the same.

For the DZ10 model, the optimal tree has now 11 leaves and it is illustrated in Fig.12. As discussed previously, the tree can be easily translated into a small numerical code using a simple structure.

In Fig.13, we illustrate the importance of the features used to build such a tree. It is interesting to observe that the most important feature is the charge dependence and the isovector dependence of the model (). By comparing with Fig.10, we observe that the relative importance of the features strongly depends on the model. In particular,

FIG. 12: Decision tree for DZ10 model given in Eq.21. For the sake of clarity and readability, the impurity (MSE) was omitted.

four features of nine turned out not to be relevant during the optimisation of the tree. We could further simplify the tree by reducing the features used or exploring new ones. This investigation goes beyond the scope of the present guide since we are only interested in illustrating how the algorithm works.

We implement such a tree within a simple Fortran code. See AppendixA for details. With such a code, we calculate the new binding energies as

The global RMS drops to = 0.471 MeV (= 0.569 MeV). The improvement on the binding energies is not as good as the one obtained [31] using a more complex neural network, but the model we produced is far simpler. It is worth noticing that the final model given in Eq.23 is fully specified by 43 parameters (the 10 original parameters from the DZ10 models, the 7 pairs describing the variable and the value to split on, the 9 values of the response on each leaf). This number is comparable with most nuclear mass models [36–39], which perform similarly to Eq.23.

As done previously for the LD model, we represent in Fig.14 the splits of the tree given in Fig.12. We see that the most important cuts take place along Z. This was also the most important feature of the model as shown in Fig.13. Interestingly, using the decision tree we have identified a large area in the residuals corresponding to the medium-heavy neutron-rich nuclei for which the correction is very small. On the contrary, the same mass range, but on the proton-rich side, requires a much more significant energy correction. This may be a symptom of poor treatment of the isovector channel in the model.

In Fig.15, we represent the comparison between the original residuals obtained with DZ10 model and the improved one using the decision tree. The histogram has been normalised. We see that the new residuals are now more clustered around the mean value, although we see that there are still some heavy tails that we have not been able to eliminate. We have checked the normality of the residual using the standard Kolmogorov Smirnov test [40] and we can say that the residuals are not normally distributed with a 95% confidence, thus showing there is still some signal left in the data that we have not been able to grasp.

We conclude this section by summarising the impact of decision trees on the residuals of the various mass models

FIG. 13: Relative importance of the features (reduction in impurity normalised) in the DZ10 model. As before, the features (defined as in the caption of Fig.10) were not used in the model, has zero importance in the model and can thus be discarded.

FIG. 14: (Colours online) Graphical representation of the splits done by the decision tree illustrated in Fig.13 on the Segr´e chart of nuclei. Corrections to the DZ10 model as a function of N, Z.

and different features used in the calculations. The results are reported in Tab.II. We observe that using feature engineering, we have been able to reduce the RMS of the LD model by 30%. Adopting the same features for the DZ10 model, we have improved the global RMS by 18%.

It is worth noting that the numbers given in Tab.II are strictly dependent on the features we used to build the trees. Different choices would lead to different numbers.

V. CONCLUSION

In this guide, we have illustrated a well-known decision tree algorithm by providing very simple and intuitive examples. We have also shown the importance of analysing the data to improve the performances of the method.

We have applied the decision tree to the case of two well known nuclear mass models: liquid drop and Duflo-Zuker. In both cases, using a small number of leaves (9 and 11 respectively), we have been able to improve the global RMS

FIG. 15: (Colours online) Normalised histograms comparison between residues with and without the tree correction.

TABLE II: Here, is the original model RMS, is the RMS once the corrections are added the RMS on unseen data (with corrections).

of the models by 29.5% and 17.6% respectively. More consistent improvements have been obtained in the literature using neural networks [31, 41, 42], but using a larger set of adjustable parameters.

We have also illustrated how to represent graphically the decision tree to better highlight the regions of the splits: this allows us to identify possible patterns in the data-set and eventually use them to improve the original model. By analysing the importance of the features, it is then possible to identify possible missing structure in the model

Finally, we have also illustrated how to translate the decision tree into a simple numerical code that could be easily added to existing ones to calculate nuclear masses.

Acknowledgments

We thank M. Shelley for helping us with this guide. This work has been supported by STFC Grant No. ST/P003885/1.

For completeness, we provide here a possible translation of the decision tree into a Fortran code.

Notice that a decision tree is formed by a simple sequence of conditional statements and the example given here can be easily ported to any other used computer language.

[1] J. Friedman, T. Hastie, and R. Tibshirani, The Elements of Statistical Learning: Data Mining, Inference, and Prediction (Springer-Verlag New York Inc., 2009), ISBN 978-0387848570.

[2] P. Mehta, M. Bukov, C.-H. Wang, A. G. Day, C. Richardson, C. K. Fisher, and D. J. Schwab, Physics reports (2019).

[3] D. G. Ireland and W. Nazarewicz, Journal of Physics G: Nuclear and Particle Physics 42, 030301 (2015).

[4] A. Pastore, Journal of Physics G: Nuclear and Particle Physics 46, 052001 (2019).

[5] B. P. Roe, H.-J. Yang, J. Zhu, Y. Liu, I. Stancu, and G. McGregor, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 543, 577 (2005), ISSN 0168-9002.

[6] H.-J. Yang, B. P. Roe, and J. Zhu, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 555, 370 (2005).

[7] S. Bailey, in New analytical techniques for the investigation of alpha clustering in nuclei (2017).

[8] L. Breiman, Statistical Science 16, 199 (2001), URL http://www.jstor.org/stable/2676681.

[9] G. E. Box, Journal of the American Statistical Association 71, 791 (1976).

[10] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, et al., Journal of Machine Learning Research 12, 2825 (2011).

[11] M. T. Ribeiro, S. Singh, and C. Guestrin, in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, August 13-17, 2016 (2016), pp. 1135–1144.

[12] O. Boz, in Proceedings of the Eighth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (Association for Computing Machinery, New York, NY, USA, 2002), KDD 02, p. 456461, ISBN 158113567X, URL https: //doi.org/10.1145/775047.775113.

[13] L. Breiman, Machine Learning 45, 5 (2001), ISSN 1573-0565.

[14] T. Chen and C. Guestrin, in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (Association for Computing Machinery, New York, NY, USA, 2016), KDD 16, p. 785794, ISBN 9781450342322, URL https://doi.org/10.1145/2939672.2939785.

[15] P. A. Lachenbruch and M. R. Mickey, Technometrics 10, 1 (1968), ISSN 00401706.

[16] C. Maderna and P. M. Soardi, Lezioni di analisi matematica (CLUED, 1985).

[17] L. Breiman, J. H. Friedman, R. A. Olshen, and C. J. Stone, Classification & Regression Trees (Wadsworth Co., 1983), ISBN 978-0534980542.

[18] G. Van Rossum and F. L. Drake Jr, Python tutorial (Centrum voor Wiskunde en Informatica Amsterdam, The Netherlands, 1995).

[19] J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).

[20] S. Milborrow, (2019), r package version 3.0.8, URL https://CRAN.R-project.org/package=rpart.plot.

[21] T. Therneau and B. Atkinson, rpart: Recursive Partitioning and Regression Trees (2019), r package version 4.1-15, URL https://CRAN.R-project.org/package=rpart.

[22] G. Louppe, arXiv 1407.7502 (2014).

[23] B. Efron, Ann. Statist. 7, 1 (1979).

[24] L. Breiman, Machine Learning 24, 123 (1996).

[25] O. Okun and H. Priisalu, in Random Forest for Gene Expression Based Cancer Classification: Overlooked Issues (2007), vol. 4478, pp. 483–490.

[26] M. Wang, G. Audi, F. Kondev, W. Huang, S. Naimi, and X. Xu, Chinese Physics C 41, 030003 (2017).

[27] A. Sobiczewski and Y. A. Litvinov, Physical Review C 89, 024311 (2014).

[28] K. S. Krane, D. Halliday, et al., Introductory nuclear physics (Wiley, 1987).

[29] J. Duflo and A. Zuker, Physical Review C 52, R23 (1995).

[30] J. Mendoza-Temis, J. G. Hirsch, and A. P. Zuker, Nuclear Physics A 843, 14 (2010), ISSN 0375-9474.

[31] A. Pastore, D. Neill, H. Powell, K. Medler, and C. Barton, Physical Review C 101, 035804 (2020).

[32] C. Qi, Journal of Physics G: Nuclear and Particle Physics 42, 045104 (2015).

[33] B. Komer, J. Bergstra, and C. Eliasmith, in Hyperopt-Sklearn: Automatic Hyperparameter Configuration for Scikit-Learn (2014), pp. 32–37.

[34] J. Heaton, arXiv 1701.07852 (2017).

[35] R. Utama, J. Piekarewicz, and H. Prosper, Physical Review C 93, 014311 (2016).

[36] S. Liran and N. Zeldes, Atomic Data and Nuclear Data Tables 17, 431 (1976).

[37] P. M¨oller, W. Myers, W. Swiatecki, and J. Treiner, Atomic Data and Nuclear Data Tables 39, 225 (1988).

[38] S. Goriely, N. Chamel, and J. Pearson, Physical Review Letters 102, 152503 (2009).

[39] N. Wang and M. Liu, Physical Review C 84, 051303 (2011).

[40] R. J. Barlow, Statistics: a guide to the use of statistical methods in the physical sciences, vol. 29 (John Wiley & Sons, 1993).

[41] R. Utama and J. Piekarewicz, Physical Review C 96, 044308 (2017).

[42] L. Neufcourt, Y. Cao, W. Nazarewicz, F. Viens, et al., Physical Review C 98, 034318 (2018).

[43] We provide an HTML version of the material at the web address https://mlnp-code.github.io/mlnp/