On the impact of selected modern deep-learning techniques to the performance and celerity of classification models in an experimental high-energy physics use case

3.1.1. Overview

3.1.2. Scoring metric

The performance of a solution is measured using the Approximate Median Significance [43], as computed in equation (1):

$$\begin{equation} \mathrm{AMS} = \sqrt{2\left(\left(s+b+b_\mathrm{reg.}\right)\ln\left(1+\frac{s}{b+b_\mathrm{reg.}}\right)-s\right)}, \end{equation} \tag{ 1 }$$

in which s is the sum of weights of true positive events (signal events determined by the solution to be signal), b is the sum of weights of false positive events (background events determined by the solution to be signal), and $b_\mathrm{reg.}$ is a regularisation term that was set to 10 by the competition organisers. This term effectively ensures that a minimum number of background events are always present in the signal region, which reduces the variance of the AMS (see section 4.3.2 of reference [39]). The AMS provides a quick, analytical value which is an approximation of the expected significance one would obtain after a full hypothesis test of signal+background versus background only. The significance is a Z score corresponding to the minimum number of standard deviations, σ, a Gaussian distributed variable would have to fluctuate in order to produce the observed data, see e.g. reference [44] for a compact introduction to hypothesis testing in new-physics searches. In High Energy Physics, an observed significance of at least five sigma is commonly required to claim observation of a new particle or process; for a discussion on the history and appropriateness of this convention one may read reference [45].

The common practice for these kinds of problems is to cut events from the data in order to remove preferentially the background events whilst retaining the signal events, in order to improve the AMS of the remaining data. Either a single cut can be used on some highly discriminating feature of the data, or multiple cuts can be used over several features. The common ML approach is to adapt the former procedure by using the features of the data to learn a new single highly discriminating feature, place a threshold on it, and then only accept events which pass the threshold. The feature learnt in this approach is simply the predicted class distribution of events, in which background events will be clustered towards zero, and signal events towards one. The AMS can then be optimised by only accepting events with a class prediction greater than some value.

In a full HEP analysis, the sensitivity is usually improved further by fitting histograms of the data along one or two discriminating features, and computing the significance using multiple bins of the data. This can either involve cutting on the ML prediction and then fitting to a different feature, or not cutting on the prediction and instead directly fitting to the histogram of the prediction. Additionally, the data may be split into analysis-specific sub-categories to further increase sensitivity, due to varying background contributions in each category. Examples of such sub-categories for this search might be: channel-wise 'tau+electron' and 'tau+muon'; and jet-wise 'zero jets', 'one jet', 'two or more jets'. However, for simplicity the challenge here requires only a single cut on the ML prediction.

The threshold cut can easily be optimised by scanning across the range of possible cuts (either at fixed steps, or checking at each data point) and picking the value which maximises the AMS. This is likely, however, to be optimal only for the dataset on which it is optimised, and performance can be expected to drop when applying the cut to unseen data. This is reflected in the challenge by requiring the solutions to predict on test data for which the class labels are not provided. It is important then that one is more careful when choosing a cut, in order to generalise better to unseen data. The approach adopted here is to consider the top 10% of events as ranked by their corresponding AMS, and compute the arithmetic mean of their class predictions and use this as the cut. This reduces the influence of statistical fluctuations in the AMS, resulting in a more generalising cut.

It is our understanding that the winning solution to the Higgs ML Challenge, with an AMS of 3.80 581, has yet to be outperformed under challenge conditions and so will be used as a comparison to the solution developed in this article. Whilst reference [46] claims to achieve an AMS of 3.979 using a method of converting the tabular data to images, we find that this result is likely a mistake in the authors' handling of the testing data which led them to compute the AMS on double the integrated luminosity, thereby increasing their score by a factor of approximately $\sqrt{2}$. We assert that their method actually achieves an AMS of around 2.8 and document our study in Appendix B.

The ATLAS detector [41] is cylindrical in shape, displaying rotational and forwards-backwards symmetry in the transverse and longitudinal planes, respectively, with respect to the beam-axis. It is worth noting that such a configuration is common amongst other particle detectors which are also designed to be suitable for searches for, and studies of, the Higgs boson, such as the CMS detector.

Two coordinate systems are normally used in particle physics: the Cartesian system of $\left(x,y,z\right)$, which can be used to express geometric position and components of particle momenta; and cylindrical coordinates $\left(r,\eta,\phi\right)$, again used for geometric position and components of particle momenta (in which case radius r is the particles' transverse momenta PT).

At ATLAS the Cartesian system is defined as shown in figure 1 such that the origin is the nominal collision-point of the experiment, the z-axis points along the beam-axis, the x axis points towards the centre of the LHC ring, and the y axis points perpendicularly upwards.

Zoom In Zoom Out Reset image size

The cylindrical coordinates for momenta are defined in terms of transverse momentum PT, azimuthal angle φ, and pseudorapidity η, where for a particle with momentum p:

$$\begin{equation} \eta = \tanh^{-1}\left(\frac{p_z}{\left|\bar{p}\right|} \right). \end{equation} \tag{ 2 }$$

In this article the ΔR separation between two objects refers to the Euclidean distance between two points in (φ, η) space:

$$\begin{equation} \Delta R = \sqrt{\Delta\phi^2+\Delta\eta^2}, \end{equation} \tag{ 3 }$$

where Δφ and Δη are the angular separations between the points in terms of azimuthal angle and pseudorapidity, remembering that φ is cyclical within [−π, π).

Although the full dataset created for the challenge has now been made public [47], in order to provide an accurate comparison to the previous scores only the subset that was used during the challenge is used here, both for training and testing. The training and testing sets consist of 250 000 and 550 000 events, respectively. Both contain a mixture of both classes: signal ($h\rightarrow\tau^{+}\tau^{-}$); and background ($t\bar{t}, Z\rightarrow\tau^{+}\tau^{-}$, and W boson decay). The full dataset is only slightly larger, containing an extra 18 238 events (an increase of 2.3%), and so it is unlikely to offer much in the way of benefits to training or evaluation. Use of the full dataset was therefore not thought to be worth the cost of being unable to compare to the baseline results on a like-for-like basis.

Each event is characterised by two sets of training features: primary - lower level information such as tau PT; and derived - higher-level information calculated via (non)-linear combinations of the low-level features or from hypothesis-based fitting procedures. Tables 1 and 2 list and describe the features, and further details are available in reference [29]. The coordinate system of the data is originally provided in terms of $\left(\textsc{PT}, \eta, \phi\right)$, however initial tests and past experience dictates that NN-based models perform better when vectors are mapped into a Cartesian coordinate system, due to the cyclical nature of the azimuthal angle, and the non-linear nature of the pseudorapidity, η.

Table 1. Glossary of high-level (derived) features used for training, along with their grouping.

Feature name	Description	Grouping
DER_mass_MMC	Mass of the Higgs boson estimated by a hypothesis based fitting technique	Mass, Higgs
DER_mass_transverse_met_lep	Transverse mass of the lepton and $\stackrel{\rightarrow}{P}_{T}^{miss}$	Mass Higgs
DER_mass_vis	Invariant mass of the lepton and the tau	Mass, Higgs
DER_pt_h	Transverse momenta of the vector sum of the lepton, tau, and $\stackrel{\rightarrow}{P}_{T}^{miss}$	3-momenta, Higgs
DER_deltaeta_jet_jet	Absolute difference in pseudorapidity of the leading and subleading jets (undefined for less than two jets)	Angular Jet
DER_mass_jet_jet	Invariant mass of the leading and subleading jets (undefined for less than two jets)	Mass, Jet
DER_prodeta_jet_jet	Product of the pseudorapidities of the leading and subleading jets (undefined for less than two jets)	3-momenta, Jet
DER_deltar_tau_lep	Separation in η − φ space of the lepton and the tau	Angular, Final-state
DER_pt_tot	Transverse momentum of the vector sum of the transverse momenta of the lepton, tau, the leading and subleading jets (if present), and $\stackrel{\rightarrow}{P}_{T}^{miss}$	Sum, Final-state
DER_sum_pt	Sum of the transverse momenta of the lepton, tau, and all jets	Sum, global event
DER_pt_ratio_lep_tau	Transverse momenta of the lepton divided by the transverse momenta of the tau	3-momenta, Final-state
DER_met_phi_centrality	Centrality of the azimuthal angle of $\stackrel{\rightarrow}{P}_{T}^{miss}$ relative to the lepton and the tau	Angular, Final-state
DER_lep_eta_centrality	Centrality of the pseudorapidity of the lepton relative to the leading and subleading jets (undefined for less than two jets)	Angular, Jet

The dataset also contains information about the hardest two jets in each event. For cases in which there were not enough jets reconstructed, the features are set to default values of -999. In order to prevent these values from having an adverse effect on the normalisation and standardisation transformations that will later be applied, these values were replaced with NaN (Not a Number - a special value that is ignored by the pre-processing methods), which prevents them influencing the transformations during fitting, and being altered by the transformations during application. Additionally, one of the features, DER_mass_MMC, which is the mass of the Higgs boson estimated by a hypothesis based fitting technique, sometimes does not converge, in which case the value is set to NaN.

Since scoring the testing data requires optimising a threshold on some data, and also to provide some way to compare architectures without relying on too much on the public score of the testing data, an 80 : 20 random, stratified split ⁵ into training and validation sets is performed on the original training data. Since it was found that the validation set had a large impact on the test scores, due to the cut optimisation, experiments are repeated multiple times using different random seeds for the splitting.

After processing the data consist of 31 training features. The training data features are transformed to have a mean of zero, and a standard deviation of one. The exact same transformation is then applied to both the validation and testing data. Following the pre-processing stage, all NaN values are replaced with zeros. The training and validation datasets are then each split into ten folds via random stratified splitting on the event class (signal or background). The testing data is also split into ten folds, but via simple random splitting. This was done to allow for easy training of ensembles and to make the process of augmenting the data easier. Additionally, we henceforth redefine one epoch as being a complete use of one of these folds, and during training will load each fold sequentially.

As mentioned in section 3.1.1, each event in the training and validation data also contains a weight. This is normalised to be used to evaluate the AMS, but also reflects the relative importance of the particular event. It is a product of the production cross-section of the underlying process and the acceptance efficiency of the initial skimming procedure. Since the model is unlikely to achieve perfect classification, it is important that it focuses on learning to classify the higher weighted events better than other, less important events. This can be achieved by applying the weights as multiplicative factors during the evaluation of the loss:

$$\begin{equation} \mathcal{L} = \frac{1}{N}\sum^N_{n = 1}w_n\times \epsilon\left(y_n,\hat{y}_n\right), \end{equation} \tag{ 4 }$$

where w_n is the weight associated with event n in a mini-batch of N events. This method of loss weighting is also used to account for class imbalances in the data by normalising the event weights in training data such that the sum of weights for each class are equal:

$$\begin{equation} \sum w_{\mathrm{signal}} = \sum w_{\mathrm{background}} = 1. \end{equation} \tag{ 5 }$$

This has the advantage of retaining the relative importance of different events within each class, whilst balancing the importance of each class overall.

The HiggsML dataset is used due to it being a publicly accessible dataset which is representative of a typical HEP problem in both size and scoring metric, and already has strong baseline solutions from the 2014 Kaggle challenge.

The public release of the dataset now includes the truth labels of the original Kaggle test set meaning that solutions can be fully scored without requiring submissions to Kaggle. This also means that the test set can be further subdivided and mean scores computed allowing both a score and an uncertainty to be reported.

The preceding version of this investigation, available in reference [30], attempted to reproduce the challenge conditions by only evaluating the public score during development and then checking the private score at the end. Additionally, that investigation only considered single trainings with a fixed seed for the train/validation splitting. The investigation presented here instead attempts to evaluate the general effects that new techniques have on model performance, and so will report the mean scores from multiple trainings, each with unique random seeds for train/validation splitting. Additionally, the mean public score over ten stratified folds will be computed.

Whilst this metric would be unavailable in true challenge conditions, the final score is a product of both model optimisation and cut optimisation. Since the method to pick the cut is fixed (as per 3.1.2), the problem of cut optimisation reduces to cut compatibility. Even so, the validation AMS alone is not representative of the true performance of the method; the score on some withheld data at the chosen cut is required. Data limitations mean that it is necessary to rely on the public scores to provide some indication of true model performance. However in order to avoid overfitting to the public score, the mean score will be used, and this will be averaged over repeated trainings.

This all is to say that, although the private scores will be kept blind until the end, this investigation does not fully reproduce challenge conditions and instead aims to demonstrate the general effect of different architectures and training choices. Additionally, the solutions developed will consider the train and inference time requirements, not just performance, and hopes to produce solutions and approaches which are appropriate for use in other HEP applications.

Several different metrics will be reported during the investigation and are averaged over six repeated trainings, each of which uses a unique train/validation splitting. The metrics are:

• The Mean Maximal Validation AMS (MMVA); the maximum AMS achievable on the validation data. MMVA can be thought of as representing an optimistic upper-limit of solution performance.
• The Mean Validation AMS at Cut (MVAC); the AMS on the validation data at the cut chosen according to the technique described in section 3.1.2. MVAC offers a more general measure of performance, but assumes that the cut used generalises well to unseen data.
• The Mean Averaged Public AMS (MAPA); the value of the AMS on the public part of the testing data at the chosen cut averaged over ten sub-splittings. MAPA represents a more realistic indication of performance since it accounts for cut generalisation.

MAPA will be the main metric for determining improvements, with MMVA and MVAC available to provide further indications of performance. Uncertainties computed for these metrics only account for statistical fluctuations in performance as systematic uncertainties due to theoretical modelling and experimental precision are not accounted for in the challenge.

Additionally, the time taken to train and apply each setup will be reported as the fractional increase in time compared to the baseline setup (see section 4.3), and this will be averaged over the six hardware setups used. A more complete set of timings will be reported towards the end of this article.

Since this paper will report the timings of the solutions, we introduce here the machinery used. There are five machines, two of which have the same configuration, and one of which is used in two different hardware configurations (CPU or GPU). Each solution is trained once per hardware configuration, meaning that six trainings are performed. Whilst RAM and VRAM capacity will be stated, we emphasise that the memory requirements of the solutions are minimal and less than 1 GB.

4.1.0.1. Desktop

4.1.0.2. MacBook Pro (2018, 13)

4.1.0.3. Virtual Machine (m2.medium)

4.1.0.4. 2 X Virtual Machine (m2.large)

The definitions of the features available, and selection of the training features to use. All training features are found to be useful.

As discussed in section 3.4, the dataset consists of both low-level and high-level features: 18 low-level features consisting of 3-momenta $(\textsc{PT},\eta,\phi)$ ⁶ , missing transverse momentum ⁷ , and sum of transverse momenta of all jets; and 13 physics-inspired high-level features such as invariant masses ⁸ and angles between final-state objects. Tables 1 and 2 describe the features used along with the associated name that will be used throughout this article. Figure 2 illustrates the density distributions of two example features, one high-level and the other low-level. It should be borne in mind that 3-momenta of final-states are the result of complex reconstruction techniques and so are not strictly 'low-level' features, these would be the tracking hits and energy deposits in the particle detector, however they are the lowest-level of information provided in the dataset.

Zoom In Zoom Out Reset image size

The mass of the Higgs boson candidate is expected to be a particularly useful feature: the τ+τ− final-state candidates should form an invariant masses around 125 GeV for signal, 91 GeV the $Z\rightarrow\tau{+}\tau{-}$ background, and a broad range of masses for the W-decay (with perhaps a peak around 80 GeV) and $t\bar{t}$ backgrounds. The tau decays, however, involve neutrinos, which cannot be detected at the LHC detectors leading to a loss of kinematical information. DER_mass_vis is an estimation of the Higgs mass using the visible decay products, and underestimates the Higgs boson mass. DER_mass_MMC instead attempts to account for the missing energy and produce a more accurate estimate for the Higgs mass by performing a kinematic fit using the visible decay products and the missing transverse momentum under the hypothesis of a di-tau decay [48].

Our hope is that the networks will be able to exploit the low-level information, but without this exploitation the low-level features risk being removed during feature selection. We are therefore more interested in seeing whether all of the high-level features are necessary. Testing will use a simplified type classification model called a Random Forest [49]. We use these as they are quick and easy to train, provide decent performance, and do not require the data to be transformed in any particular fashion prior to training.

Quantification of a feature's importance is done using its permutation importance. This aims to quantify how much a given feature influences the model's output. Features with a higher importance will have a greater influence overall that those of lesser importance. The permutation importance (PI) is computed by first training the model as normal and evaluating its performance on the training set. Next, a copy of the training data is made and one of the features is randomly shuffled datapoint-wise, thereby destroying any information carried by that feature for each datapoint, whilst keeping its values at the same scale as seen during training. The performance of the model on the shuffled data is computed and the difference in performance compared to the original performance is the permutation importance of the feature that was shuffled. If shuffling the feature caused a large change in performance, then that feature was important to the model. If the change was low, then the feature was less important. By then dividing the PI by the original performance we arrive at the fractional permutation importance.

This process of copying and shuffling is repeated for every feature. By retraining new models, an average PI can be computed for each feature. Although the dataset is the same each time, the feature subsampling at each split acts as a source of stochasticity. Here we compute PIs using an average of five retrainings. A threshold on the PI can then be used to select features by minimum importance. We use a threshold of 0.005 as to identify features as 'important'. From figure 3, we can see that all features except the four nearest the top of the plot, i.e. those associated with jet pairs, appear to be important.

Zoom In Zoom Out Reset image size

The lack of importance of jet-pair features is to be expected: since the majority of events do not contain two or more jets, the values for these features are often the default value of zero. When the feature is permuted, most of the zeros will be replaced with another zero, and the model response remains mostly unaffected. Effectively the di-jet features are only important for a subset of events, but their importance is smeared out across all events. In actual application, the jet number will be encoded in a way such that the model will be able to more easily target these di-jet events and make use of their features. The di-jet features will therefore be kept. As additional motivation, training a new Random Forest without these features showed a drop in performance on unseen data.

In order to check for monotonic relationships between features we then compute the Spearman's rank-order correlation coefficients [50] for every pair of high-level features and do not find any pairs that are fully (anti-)correlated with one another. As a final check we examine the mutual dependencies of the features, i.e. how easy it is to predict the value of a feature using the other features as inputs to a Random Forest regressor. As shown in figure 4, none of the non-di-jet high-level features are completely predictable using the other features (dependencies are never equal to one). We can therefore expect that each high-level feature is bringing at least some unique information not carried by the others.

Zoom In Zoom Out Reset image size

The basic architecture for later comparisons, and also a demonstration of the benefits of ensembling models.

Training and inference of all solutions developed in this article are performed in LUMIN [51], a wrapper library for PYTORCH [52], which so far has been primarily developed by the author of this article.

4.3.1. Architecture

4.3.2. Performance

4.3.2.1. Single model

A single model is trained on nine out of the ten training folds, with the tenth fold acting as validation data for loss monitoring. Once finished, the model is applied to both the validation and testing data, and the required metrics recorded. This process is run once per hardware configuration, as defined in section 4.1 (six times in total), each using a unique random seed for the train/validation splitting as described in section 3.4. Performance metrics for the single model are reported in table 3.

Table 3. Validation-performance comparison between using the single model and weighted ensemble setups. Higher values of MMVA, MVAC, and MAPA indicate better performance, and higher values of fractional time-increase indicate that the solution takes longer to train and/or apply. The best values for each metric are shown in bold, and the setup chosen is also indicated in bold.

				Fractional time-increase
Setup	MMVA	MVAC	MAPA	Training	Inference
Single model	3.43 ± 0.04	3.38 ± 0.04	3.44 ± 0.02	–	–
Ensemble	$\mathbf{3.70\pm0.05}$	$\mathbf{3.64\pm0.04}$	$\mathbf{3.664\pm0.007}$	10.4 ± 0.6	5.9 ± 0.9

4.3.2.2. Ensemble

An ensemble of ten models is trained, again using nine out of the ten training folds, with the tenth fold acting as validation data. Each model uses a different fold of the training data for validation. The resulting models are ensembled by weighting their predictions proportionally to the inverse of their associated loss on their validation folds. The prediction of the ensemble is then:

$$\begin{equation} p_{\mathrm{ensemble}}(x) = \sum_{i}^{10} w_ip_i(x), \end{equation} \tag{ 6 }$$

where $p_{\mathrm{ensemble}}(x)$ is the prediction of the ensemble for datapoint x, p_i (x) is the prediction of model i for datapoint x, and w_i in the importance weight of model i. Importance weights are normalised such that they sum to one. This means that models that reached a lower validation loss during training have a higher importance weight and consequent a greater influence over the ensemble prediction than those with a higher loss.

Ensembling is seen to provide improvements in all metrics compared to the single model setup, with expected increases in train and inference time. Note that the inference-time increase appears to be sub-linear due to the fact that the data are kept in memory between model evaluations. Based on the improved performance, we proceed to adopt the ensembled version as the baseline model for later comparisons.

The testing of an efficient way to encode categorical information, which results in mild improvements.

In section 4.3 we considered all features in the data to be continuous variables, however the number of jets in an event (PRI_jet_num) could instead be thought of as a categorical feature. Although two jets are more than one jet, this can also be treated as defining subcategories of particle interactions, e.g. a two-jet event and a one-jet event, in which case there is no longer a numerical comparison. Optimal treatment of categorical features requires encoding the different categories in a way such that the ML algorithm can easily access them individually. A common way to do this is to assign each category an integer code which begins at zero and sequentially increases, and then treat these codes as row indices of an embedding matrix. The elements in the indexed row are then used to represent the corresponding category.

Under the 1-hot encoding scheme, the embedding matrix is a N ×  N identity matrix (where N is the number of categories in a feature; its cardinality), i.e. each category is represented by a vector of length N in which only one value is non-zero (one element is hot). Whilst this approach provides unique access to each category, the number of input variables is now N. For high cardinality features, or many low-cardinality features, this can quickly lead to a huge increase in the number of model parameters which must be learnt.

Reference [59] instead suggests for neural networks to use an N ×  M weight matrix for the embedding, where $M \lt N$ is a hyper-parameter which must be chosen. Each category is now represented by a real-valued vector of length M, i.e. in a more compact fashion than 1-hot encoding would provide. The values of the weight matrix are then learnt via backpropagation during training. This approach of categorical entity embedding works because, for example, it is likely that it is more important to know whether an event contains zero, one, several, or many jets, rather than exactly how many jets, but rather than manually compacting the categories (and having to define a priori 'several' and 'many'), a compact, information rich, and task specific representation can be learnt. The additional hyper-parameter M can be estimated using domain knowledge or by using some rule-of-thumb: the original paper suggests to use M = N − 1, and FastAI suggests to use half the feature cardinality or 50, which ever is smaller ($M = \textrm{min}\left(50, \left(N+1\right)//2\right)$ [60])

For embedding PRI_jet_num, M is chosen to be three and initial weights for the 4 ×  3 matrix are sampled from a unit-Gaussian distribution. As shown in table 4, this is found to provide metric improvements across the board, for only minor increases in train and inference time. Figure 5 shows an example of one of the learned embeddings. Please be aware that table 4 includes all improvements that are discussed in this article, and so contains solutions which may not have been introduced at time of reading.

Zoom In Zoom Out Reset image size

Whilst not required for this particular study, as discussed in section 3.1.2, HEP analyses potentially categorise data in order to further improve the sensitivity of their search. These inference categories are search channel-dependent, such as di-tau decay-channel, as well as categorisation by jet-multiplicity, but category flags can potentially be encoded as categorical features and embedded. This then allows the classifier to better parameterise its responses to each inference category, whilst benefitting from a larger training sample (as opposed to training separate classifiers for each category) and learning information common to all categories.

Discussion and application of HEP-specific data augmentation, which result in large improvements.

4.5.1. Overview of data symmetries

At the LHC [32], beams of protons are collided head-on with approximately zero transverse momentum. Because of this the resulting final states can be expected to be produced isotropically in both the transverse plane (x, y) and the longitudinal axis (z). Particle detectors such as CMS [61] and ATLAS [41] are built to account for these isotropies by being symmetric in both azimuthal angle (φ) and pseudorapidity (η). All this is to say that the class of physical process that gave rise to the collision event is only related to the relative separation between the final states, and not the absolute position of the event.

Since the data used in this problem are simulated for such a collider and detector combination, the class is invariant to flips in the transverse or longitudinal planes, and to rotations in the azimuthal angle, as illustrated in figure 6. Note, it is important to reconsider the set of class-preserving transformations if intending to use the following techniques to analyse data for a different experimental configuration, such as: colliding beams of different particles, using a fixed-target experiment, or an asymmetrical detector.

Zoom In Zoom Out Reset image size

4.5.2. Data fixing

4.5.3. Data augmentation

4.5.4. Comparison

Table 5 compares the two approaches in terms of the optimisation metrics. Data augmentation produces more favourable scores, but at the expense of a large increase in inference time. Nevertheless, the absolute inference time is still acceptable (15 s on a GPU (Nvidia GeForce GTX 1080 Ti) and 2 min on the slowest CPU (Intel Xenon Skylake CPU, 2 virtual single-thread cores, clock-speed 2.2 GHz)), so the augmentation approach was chosen going forwards. Table 4 compares the improvements against the baseline solution.

Table 5. Comparison of data symmetry exploits in terms of the optimisation metrics and time impact. The best values for each metric are shown in bold, and the setup chosen is also indicated in bold.

				Fractional time-increase
Setup	MMVA	MVAC	MAPA	Training	Inference
Fixing	3.90 ± 0.04	3.83 ± 0.05	3.76 ± 0.01	–	–
Augmentation	$\mathbf{3.96\pm0.04}$	$\mathbf{3.89\pm0.05}$	$\mathbf{3.79\pm0.01}$	1.11 ± 0.05	13 ± 5

Testing of learning-rate annealing schedules. Whilst this results in moderate improvements it eventually becomes obsolete by the schedule in section 4.9.

According to reference[63], the learning rate (LR) is one of the most important parameters to set when training a neural network. Whilst initial, optimal values can be found using a LR Range Test [58], it is unlikely that this will remain optimal throughout training. A common approach is to decrease the LR in steps [64] during training (either at preset points or when the loss reaches a plateau).

Reference [65] instead proposes a cosine-annealing approach in which the LR is decreased as a function of mini-batch iterations according to half a cycle of a cosine function. This means that initial training is performed at high LRs, where the model can quickly move across the loss surface towards a minima, and the final training is spent at lower LRs allowing the model to converge to a minima, with a smooth, approximately linear transition between these states. Additionally, by repeating this cosine decay, the paper suggests that the rate of convergence can be further improved, a process described as warm restarts.

Reference [66] suggests these sudden changes in the LR allow the model to explore multiple minima. This exploration of minima could in turn lead to the discovery of wider, less steep minima, which should generalise better to unseen data since changes in the centre of the minima will have less of an effect on the loss than if the minima were steeply sided.

Figure 7 illustrates the LR schedule and an example of the validation-loss evolution.

• (a)  
  We begin at an initial LR, γ₀, of 2× 10⁻³ and an initial cycle length, T, of one complete use of the training folds.
• (b)  
  During the cycle, the LR is decreased according to $\gamma_t = 0.5\,\gamma_0\times \left(1+\cos\left(\pi\,t/T\right)\right)$, where t is incremented after each minibatch update
• (c)  
  Once the LR would become zero (t = T), t is reset to zero, resulting in a sudden increase in the LR (the warm restart). Additionally, the cycle length T is doubled allowing the model to eventually spend longer and longer time at lower LRs, improving convergence.

Zoom In Zoom Out Reset image size

The warm restarts are visible in the validation loss, appearing as small increases in the loss. Initially, the early stopping criterion was modified such that training stopped once the model went one entire cycle without an improvement in the validation loss, however since testing showed models always reached this point in the same annealing cycle, the number of training epochs was instead fixed to avoid unnecessary training time.

The results (table 4 Row '+SGDR') show a minor improvement in both MVAC and MAPA, and a decrease in MMVA, indicating a less overly optimistic solution with better generalisation.

The inclusion of a newer activation function, which provides slight improvements.

Whilst ReLU is a good default choice for an activation function it does exhibit several problems, which more recent functions attempt to address:

• The Parametrised ReLU (PReLU) [53] is similar to ReLU, but can feature a constant, non-zero gradient for negative inputs, the coefficient of which is learnt via backpropagation during training;
• The Scaled Exponential Linear Unit (SELU) [67] uses careful derived scaling coefficients to keep signals in the network approximately unit Gaussian, without the need for manual transformation via batch normalisation [68]. It is recommended to use the lecun initialisation scheme [69] for neurons using the SELU activation function;
• The Swish function [70] was found via a reinforcement leaning based search, and features the interesting characteristic of having a region of negative gradient, allowing outputs to decrease as the input increases: $\mathrm{Swish}\left(x\right) = x\cdot\mathrm{sigmoid}\left(\beta x\right)$, where β can either be kept constant or be learnt during training. The recommended weight initialisation scheme for Swish is the same as the one for ReLU, i.e. He.

A comparison of these functions is shown in figure 8.

Zoom In Zoom Out Reset image size

Previous testing, documented in reference [30], found that the Swish activation function provided superior performance, and so its improvements were again tested here. We use the Swish-1 version, in which β = 1 and is not changed during training. As shown in Row '+Swish' of table 4, its inclusion provides slight improvements to MAPA, and so its use was accepted. Note that the small time increase is due to the exponentials in Swish being more complicated to compute than the max operation of ReLU.

A discussion and experimental test of several more complicated methods of building an ensemble of models. Additionally a new implementation of Stochastic Weight Averaging is introduced which does not rely on a predefined starting point. Whilst this provides slight improvements in training time it is eventually made obsolete by the schedule in section 4.9.

As mentioned in section 4.6, reference [66] suggests that reference [65]'s process of restarting the LR schedule allows the training to discover multiple minima in the loss surface. The paper builds on this by introducing a method in which an ensemble of networks results from a single training cycle by saving copies of the networks just before restarting the LR schedule, i.e. when the model is in the minima. This process of Snapshot Ensembling (SSE) is further refined into Fast Geometric Ensembling (FGE) in reference [71], by forcing the model evolution along curves of constant loss between connected minima (simultaneously and independently discovered by reference [72]).

The problem with these methods is that whilst they allow one to train an ensemble in a much reduced time, one still has to run the data through each model individually, so the application-time is still increased. Reference [73] instead finds an approach which leads to an approximation of FGE in a single model. This is done by averaging the models in weight space, rather than in model space. The general method involves training a model as normal, until the model begins to converge. At this point the model continues to train as normal, but after each epoch a running average of the model's parameters is updated.

In the training, one has to decide on when to begin collecting the averages: too early, and the model average is spoiled by ill-performing weights; too late, and the model does not explore the loss surface enough to allow SWA to be of use. Additionally, SWA may be combined with a cyclical learning rate, in which case weight averaging should take place at the end of each cycle.

4.8.1. Implementation

4.8.2. Comparison

Comparing the approaches to the current best solution in table 6, we can see that no approach is able to improve the current MAPA, however SWA with a constant LR achieves the same score in a shorter train time. From an example plot of the validation loss over time (figure 9) we can see that SWA not only demonstrates better performance, but it also shows a heavy suppression of loss fluctuations and converges to a loss plateau, where as the non-averaged model eventually over-fits (signified by an increasing loss in the later stages of training). The sharp drops in SWA loss are due to the old average being replaced with the newer average.

Zoom In Zoom Out Reset image size

Table 6. Comparison of the various advanced ensembling approaches. 'Current solution' refers to the solution as of section 4.7. The best values for each metric are shown in bold, and the setup chosen is also indicated in bold.

				Fractional time-increase
Setup	MMVA	MVAC	MAPA	Training	Inference
Current solution	3.96 ± 0.06	3.86 ± 0.05	$\mathbf{3.81\pm0.02}$	-	–
SSE	3.95 ± 0.05	3.85 ± 0.05	$\mathbf{3.81\pm0.02}$	0.15 ± 0.03	1.5 ± 0.3
FGE	3.95 ± 0.06	3.86 ± 0.05	3.80 ± 0.02	$\mathbf{-0.13\pm0.03}$	6 ± 1
SWA constant	3.96 ± 0.06	3.88 ± 0.05	$\mathbf{3.81\pm0.02}$	−0.09 ± 0.01	–
SWA linear cycle	$\mathbf{3.97\pm0.06}$	3.87 ± 0.04	$\mathbf{3.81\pm0.02}$	0.25 ± 0.03	–
SWA cosine	3.95 ± 0.05	$\mathbf{3.89\pm0.05}$	$\mathbf{3.81\pm0.02}$	0.22 ± 0.04	–

The application of a schedule for both learning rate and momentum, which provides a significant reduction in training time.

Reference [74] introduces the concept of super convergence, in which a specific LR schedule, 1cycle, is used to achieve convergence much quicker (between five and ten times) than traditional schedules allow. This is further discussed in reference [75]. The 1cycle policy combines both cyclical LR and cyclical momentum ⁹ (or β₁ in case of ADAM), and evolves both hyper-parameters simultaneously and continuously from their starting values to another value and then back again over the course of training in a single cycle. Hyper-parameter evolution takes place in opposite directions, with the LR initially increasing and the momentum initially decreasing, i.e. when the LR is at its maximum, the momentum is at its minimum, as illustrated in figure 10. In this way the two parameters help to balance each other, allowing the use of higher learning rates without the network becoming unstable and diverging. In the final stages of training, the LR is allowed to decrease several orders of magnitude below its initial value.

Zoom In Zoom Out Reset image size

One of the key points of reference [75] is that super-convergence can only occur if the amount of regularisation ¹⁰ present in the setup is below a certain amount, but that once super-convergence can be achieved, other forms of regularisation can then be tuned accordingly. In our current setup, we have avoided explicitly adding tuneable forms of regularisation ¹¹ for this reason. However section 4.5's choice to use data augmentation rather than data fixing equates to choosing not to reduce regularisation through data complexity. Because of this, we will be conservative in our testing of super-convergence and attempt to only halve our training time.

Implementation of 1cycle follows the suggestion of fast.ai [79], in which the halves of the cycles are half-cosine functions, rather than linear interpolations. This can be expected to offer the same advantages of the cosine annealing schedule of reference [65] and also provide a smooth transition between the directions of parameter evolution. The total length of the cycle was set to 135 folds (approximately half that required for the solutions of section 4.6 and section 4.8). The lengths of the first and second parts of the cycle were set to 45 folds and 90 folds, respectively, i.e. a ratio of 1:2, allowing a longer time for final convergence. Based on a LR range test performed at a constant β₁ of 0.9 and the knowledge that β₁ would be at a minimum when the LR is maximised, a slightly higher maximum LR of 1× 10⁻² and an initial value of 1× 10⁻⁴ were chosen.

As shown in figure 10, during the first part of the cycle, the LR is increase from its initial value (1× 10⁻⁴) to its maximum of 1× 10⁻² while β₁ is decreased from 0.95 to 0.85. In the second half, β₁ returns back up to 0.95, but the LR is allowed to tend to zero. As shown in Row '-SWA +1cycle' of table 4, although 1cycle provides the same level of MAPA performance as SWA, it is able to reach this level of performance in half the time.

An investigation into different connections between layers, which is found to provide a mild improvement in performance whilst reducing the number of free parameters in the networks.

The fully-connected architectures used so far aim to learn a better representation of the data at each layer, however this means that potentially useful information is lost after each layer if it has not been sufficiently well represented by the latest layer. Being able to do so requires that the layers have a sufficient number of free parameters to avoid losing useful information. Due to limited training data, this limits the depth of the networks which can be trained and so potentially limits the performance of the model.

Consider an example as illustrated in figure 11 for a set of inputs a, b, c, d, e, f in the case in which it is optimal to learn first a new feature based on inputs a, b, c, d, e using hidden layer 1 (e.g. the angle between the tau and the electron or muon) and then use hidden layer 2 to combine this new feature with input f (e.g. the angle multiplied by the amount of missing transverse energy). Hidden layer 1 must now learn either an identity representation of input f or a monotonic transformation of it, which requires 'budget' from both the available free parameters and the amount of training data. It would be more efficient if instead hidden layer 2 were to have direct access to input f. Consider further the case in which it is hidden layer 3 that requires access to input f, now both hidden layer 1 and 2 must learn identity representations of f.

Zoom In Zoom Out Reset image size

In their 2016 paper, 'Densely Connected Convolutional Networks' [80], Huang and Weinberger present an architecture in which, within subgroups of layers, the outputs of all previous layers are fed as inputs into all subsequent layers via channel-wise concatenation. This means that lower-level representations of the data are directly accessible by all parts of the block. This potentially allows both a more diverse range of features to be learnt, and for layers to be trained via deep supervision [81] due to more direct connections when back-propagating the loss gradient. Additionally, it means that the weights which were previously required to encode the low-level information may no longer be necessary.

Whilst the paper presents this dense connection in terms channel-wise concatenation of the outputs of convolutional filters, the same idea can be applied to fully connected networks by concatenating width-wise the output tensor of each linear layer with the that layer's input tensor; i.e. $x_{i+1} = D_i\!\left(x_i\right) \oplus x_i$, where x_i is the hidden state before the $i^{\mathrm{th}}$ linear layer (D_i ), as illustrated in figure 12.

Zoom In Zoom Out Reset image size

This concept can also be thought of as accumulating a tensor of hidden states, as illustrated in figure 13, in which each hidden layer takes as input the current state of the accumulated tensor and then concatenates its output to the tensor. In this way every subsequent hidden layer has direct access to all previous representations of the inputs, meaning that no identity representations need to be learnt and hidden layers are more easily optimised due to the more direct paths for the loss gradient during back-propagation. As discussed in reference [64], in the context or residual connections, such 'skip connections' also help mitigate the effects of over-parameterising a model; rather than forcing redundant layers to learn identity mappings, they can instead be skipped (ignored) or pushed towards a zero output. This can be expected to offer some level of leniency when choosing the number of hidden layers to used in the network.

Zoom In Zoom Out Reset image size

4.10.1. Testing

The architecture up to now has used four hidden layers, each with 100 neurons. Including the embedding and output layers, the total number of free parameters is 33 813 ¹² . For the densely connected architecture we reduce the width of linear layers to 33 (the number of continuous inputs plus the width of the categorical embedding output), and increase the depth of the network to six hidden layers. The outputs of each hidden layer except the last are concatenated with their inputs to compute the input tensor to the next hidden layer. This results in a total of 23 113 free parameters, ¹³ i.e. about two thirds of the original size. This architecture is illustrated in figure 14.

Zoom In Zoom Out Reset image size

Training of the networks used the same 1cycle schedule as presented in section 4.9. As shown in Row '+Dense' of table 4, the densely connected architecture provides mild improvements in both MVAC and MAPA, despite being significantly smaller, however due to the extra depth and concatenation operations, training and inference timing is approximately the same as the previous architecture.

Unfortunately, a mistake entered into the code during this experiment that led to the ensemble of trained models being uniformly weighted, rather than performance weighted (as per section 4.3.2.2) as intended. This mistake was only discovered during an ablation study performed at the end of the analysis after the private AMS values had been computed. In order to reflect the information seen during the course of the study, and not to draw conclusions from the private AMS, table 4 shows the performance of the uniformly-weighted ensemble, which achieves a MAPA of 3.82 ± 0.02. A retraining using performance-weighting achieves a MAPA of 3.81 ± 0.01 (comparable with the 3.81 ± 0.02 achieved by the solution of section 4.9), which might have caused us not to use dense connections: use of dense connections is slightly slower to train, but the slight reduction in variance of performance would have been seen as favourable. This mistake remained in the code throughout section 4.11 and Appendix A, meaning that the comparisons of metrics to what was seen in this section are still valid.

In this section we perform a search over the remaining hyper-parameters to finalise our solution. Interestingly, though, the existing architecture is found to work best and the model displays only a weak dependence on its hyper-parameters.

4.11.1. Overview of remaining hyper-parameters

4.11.2. Parameter scan

While some choices have allowed a reduction of the number of dimensions of the parameter search space, the remaining space still has five dimensions (depth, width, dropout rate, L₂, and growth rate), meaning a full grid-search of all possible combinations is unfeasible; instead a random search will be used. For each architecture, the following parameter-sampling rules were used:

• Depth sampled from $[2,9)\in\mathbb{Z}$
• Width sampled from $[33,101)\in\mathbb{Z}$
• Growth rate sampled from $[-0.2,1)\in\mathbb{R}$
• L₂ sampled from $\{0,10^{-2},10^{-3},10^{-4},10^{-5},10^{-6}\}\times 10^{-5}$
• Dropout rate sampled from {0, 0.05, 0.1, 0.25, 0.5}
• The total number of free parameters must be in the range $(5\times 10^{-3},1\times 10^{-5})$
• In the case of negative growth rate, the width of the last hidden layer before the output layer must be greater than one.

Just over 350 parameter sets were tested. Each test consisted of training an ensemble of three networks with the specified architecture. The MVAC and MAPA were then calculated and recorded. Figure 15 illustrates the distributions of the sampled parameters. Figure 16 shows MAPA and MVAC as a function of the total number of free parameters in the networks. From the flat linear-fit, we can see that performance is unrelated to model complexity, and that whilst the data and training methodologies are sufficient for optimising large models, such models are unnecessary to reach better performance. This is possibly due to the leniency on over-parameterised networks discussed in section 4.10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The testing results were then used to fit Gaussian process models in order to evaluate the partial dependencies of the metrics on the hyper-parameters. Figure 17 illustrates the one- and two-dimensional partial dependencies of MVAC and MAPA on the hyper-parameters. Note that the machinery used expects a minimisation problem, hence why the negatives of the metrics are shown.

Zoom In Zoom Out Reset image size

From the results, the trends in the partial dependencies of both MVAC and MAPA appear to favour minimal explicit regularisation. In terms of depth and width, MVAC shows weak dependence on width and growth rate, and slight dependence on depth, with networks with fewer than four hidden layers performing poorly compared to deeper networks. MAPA again shows weak dependence on width and growth rate, but demonstrates a bounded minimum for depth, favouring networks with around five hidden layers.

4.11.3. Further testing

Based on the scan results and partial dependencies, a number of promising parameter sets were further explored with complete trainings, and the LRs set using LR range tests. The performances of these architectures are compared to the current solution (depth 6, width 33, growth rate 0, and no dropout or weight decay) in table 7. As can be seen, none of the new architectures improve on the current architecture in terms of MAPA and although some are able to reach the same level of performance they do not provide significant reductions in terms of training or inference time.

Table 7. Comparison of the various explored architectures. 'Current' refers to the solution as of section 4.10. Architectures names are of the form: (depth) (width) (growth rate). The best values for each metric are shown in bold, and the setup chosen is also indicated in bold.

				Fractional time-increase
Setup	MMVA	MVAC	MAPA	Training	Inference
Current	3.95 ± 0.04	3.89 ± 0.05	$\mathbf{3.82\pm0.02}$	-	-
4 100 0.1	3.9 ± 0.05	3.83 ± 0.05	3.74 ± 0.01	0.23 ± 0.09	0.6 ± 0.1
4 40 0.45	3.96 ± 0.05	3.90 ± 0.05	$\mathbf{3.82\pm0.01}$	$\mathbf{-0.02\pm0.05}$	$\mathbf{-0.01\pm0.03}$
5 100 -0.1	$\mathbf{4.03\pm0.05}$	$\mathbf{3.91\pm0.05}$	3.80 ± 0.02	0.15 ± 0.07	0.6 ± 0.1
5 33 0.9	3.95 ± 0.05	3.89 ± 0.04	3.81 ± 0.02	0.19 ± 0.08	0.6 ± 0.2
8 33 -0.2	3.95 ± 0.06	3.89 ± 0.06	$\mathbf{3.82\pm0.02}$	0.4 ± 0.1	1.2 ± 0.4

Having built up a solution using several optimisation metrics and timing changes as guides, we now come to checking the performance of the solution using the remaining two metrics: The overall AMS on the public and private portions of the testing data (reminder: MAPA was only an average over subsamples of the public portion of the testing data). In figure 18 we show all metrics as a function of solutions developed, averaged over the six trainings that were run per solution. From the results, we can see that whilst MAPA consistently over-estimated actual performance, its general trend is nonetheless representative of changes in the public and private AMS.

Zoom In Zoom Out Reset image size

Also of interest for real-world application is the time taken by each solution. Plotted in figure 19 and figure 20 are the absolute (wall clock) and relative timings for each solution for the variety of hardware tested, as computed by timeit. Note that one solution was run per hardware configuration, except for the quad-core Xenon setup, which shows the mean of two different machines with the same hardware specification. The numbers in parentheses indicate the number of cores and number of threads per core, respectively, except for the Xenon machines, which are cloud-based and are allocated virtual CPUs. From these results, we can see that whilst the GPU is clearly superior, the architectures run sufficiently well on CPU to be viable for use in a particle-physics analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is also interesting to note the fractional test-time increase for the GPU, which appears to show sub-linear scaling when ensembling, but poor scaling when running inference-time augmentation. This is because LUMIN currently loads each data-fold from the hard-drive into either RAM or VRAM sequentially, and the timing includes this. Once in memory, however, all models in the ensemble are run over the fold. Since the load time for VRAM is non-negligible compared to the prediction time, the fractional increase when evaluating the ten models is much less than nine. Additionally, LUMIN currently applies data-augmentation when loading the data fold, rather than augmenting the data in-place, meaning that the load time is incurred an extra seven times when running TTA, hence the apparent poor scaling when running TTA on a GPU.

Of the solutions developed, the one from section 4.10 (densely connected network with 1cycle training scheme) offered the highest value for MAPA seen during solution development, and is also much quicker to train compared to solutions that used constant or cosine-annealed learning-rate schedules. Its performance with unweighted ensembling is public : private AMS (3.80 ± 0.01 : 3.803 ± 0.005), and with weighted ensembling (3.79 ± 0.01 : 3.806 ± 0.005), see section 6.2 for further discussion on this.

As can be seen from the architecture-search results, no explicit regularisation was needed for this problem, and performance was only slightly dependent on the widths of the hidden layers. The rule-of-thumb used to set the initial layer width for the densely connected network (width number = of inputs features (continuous + categorical embeddings)), seems to be appropriate, but further work will be necessary to see how well this rule generalises to other problems and applications.

A comparison to the solutions developed during the running of the challenge is presented in table 8. In terms of hardware, for the solutions developed here, the GPU results are for a Nvidia GeForce GTX 1080 Ti, and the CPU result is for a Intel Core i7-8559U (2018 MacBook Pro). RAM and VRAM usage for each setup was found to be less than 1 GB. Details of the other solutions are as follows:

• (a)  
  The 1^st place solution (available at https://github.com/melisgl/higgsml) primarily used a setup with a Nvidia Titan GPU and 24 GB RAM, but reference [83] also details timings for an Amazon EC2 m2.4xlarge instance (8-vCPU & 15 GB RAM);
• (b)  
  The 2^nd place solution (https://github.com/TimSalimans/HiggsML) was run on an Amazon EC2 m2.4xlarge instance (8-vCPU & 64 GB RAM) and notes the high RAM requirements of the training;
• (c)  
  And the 3^rd place solution used a 2012 quad-core laptop (https://www.kaggle.com/c/higgs-boson/discussion/10 481).

The difference in hardware makes a direct comparison of solution timing difficult, but we can try to scale the timings: The 1^st place solution used a Nvidia Titan GPU (2013, 1500 gigaflops processing power in double precision) and our solution uses a Nvidia 1080 Ti (2017, 10 609 gigaflops processing power in single precision). Accounting for the changes in processing power, the 1^st place solution should take around 100 min to train on our GPU. Our solution therefore provides an estimated speedup of 92% on GPU and 86% on CPU (comparing to GPU). Given the uncertainties involved in the metric it would be reasonable to say that the solution developed here achieves comparable performance to winning solutions of the challenge in just a fraction of the time.

Whilst the Kaggle solutions only report the private AMS scores they achieved (single numbers with no uncertainties) it is of interest to know the general performance offered by each of the solutions. The summary paper for the challenge (reference [39]) estimates the statistical significance of the solution ranking by bootstrap resampling of the scored data. The authors find that at a 95% confidence level, the 1^st -place solution does indeed outperform the rest, but that the 2^nd and 3^rd -place solutions are indistinguishable in terms of performance. Ideally such a test should also be performed to compare our solution to the 1^st -place solution, however this requires access to the scored data that was submitted during the competition, which is not publicly accessible.

Due to all possible permutations of new methods not being considered, it is plausible that some methods would show a greater or worse relative improvement if testing had been performed in a different order. This could happen for a number of reasons, such as: overlaps in improvements in methods; certain pairs of methods working well together; and the non-linear nature of the scoring metric. Table 9 addresses some these reasons by detailing results from retrainings of solutions in which only one new method was added. Whilst all new methods provide at least minor improvements in isolation, it is interesting to see that data augmentation actually provides a greater improvement than ensembling when applied individually, in comparison to figure 21 (which shows the relative contributions of each new method to the overall improvement in performance), indicating an overlap in the two methods. Also the use of dense connections on their own provides a smaller benefit than they do in the full solution, implying that they are most useful when applied in conjunction with other methods.

Zoom In Zoom Out Reset image size

In section 4.3.2.2, performance-weighted ensembling was used. The more simple form of ensembling would have been to instead consider the predictions of each model equally (uniform weighting). The choice to use weighted ensembling was made based on previous experience with similar HEP classification tasks, and was not tested at the time so as to reduce the number of times the performance metrics were computed (reducing the chance of over-fitting to the validation dataset). Having completed the study, though, it is of interest to check the impact of weighted versus unweighted ensembling. As discussed in section 4.10, the final ensembles were accidentally built using uniform weighting due to a mistake in the code, achieving a private AMS of 3.803 ± 0.005. The mistake was found when performing this ablation study of the ensemble weighting and retraining the ensembles using performance-weighted ensembling, as intended, results in a slightly higher private AMS of 3.806 ± 0.005.