Machine learning for particle physics using R – Data production pipeline – Visualising a correlation matrix as a force-directed network

Quark/gluon jet tagging for ALICE

Machine learning for particle physics using R

Andrew John Lowe

Wigner Research Centre for Physics,

Hungarian Academy of Sciences

Introduction: about this talk

• This is a walk-through of an ongoing analysis that I am performing entirely in R
• Computing power was limited to my own laptop
• Limited RAM and CPU processing power imposed contraints that drove the analysis decisions
• Similar constrains apply when performing machine learning in TMVA (ROOT) or Scikit-learn (Python), both of which are used in HEP, so the lessons learned are applicable when using these tools also
• This talk focuses on the problem of distinguishing between (light) quark- and gluon-initiated jets

The problem in a nutshell

• Beams of energetic protons and/or heavy ions collide inside our detector
• Quarks and gluons emerge and decay into collimated sprays of particles
• Algorithms cluster these decay products into jets
• For each jet, we'd like to know what initiated it
• Was it a Quark? Or a Gluon?
• Being able to accurately discriminate between quark- and gluon-initiated jets would be an extremely powerful tool in the search for new particles and new physics
• This is an archetypal classification problem that might be amenable to machine learning

Data production pipeline

• Use experiment's software (AliROOT) to process Monte Carlo simulated data that contains lots of jets
  • Insert my own C++ code with handcrafted features
• Attach ground-truth class labels ("quark"/"gluon") to data for each jet
  • There is significant class noise (mislabelled jets)
  • Different labelling schemes exist, but none are perfect because Monte Carlo simulation cannot perfectly simulate real data
• Write-out data and convert to CSV format for use in R

Class labelling

• Jets labelled using the partons in the generator event record
• Parton with highest \(p_{\mathsf{T}}\) within a \({\Delta}R\) equal to the radius parameter of the jet algorithm determines the jet label
  • This is identical to the scheme used by ATLAS\(^{*\dagger}\)
  • This labelling procedure is not unambiguous and is not strictly identical for different MC generators
  • Definitions are not theoretically robust, but studies (with MADGRAPH) have shown that for most generators, truth labelling is identical to matrix-element-based labelling for 90-95% of (isolated) jets

* Light-quark and gluon jet discrimination in \(pp\) collisions at \(\sqrt{s}\) = 7 TeV with the ATLAS detector, Eur.Phys.J. C74 (2014) 3023

\(\dagger\) Light-quark and Gluon Jets in ATLAS: Calorimeter Response, Jet Energy Scale Systematics, and Sample Characterization, ATLAS-CONF-2011-053, Mar. 2011, also ATLAS-CONF-2012-138, Sept. 2012

What about mislabelled jets?

• Impose a requirement that the jets are isolated to restrict contamination from wide-angle QCD radiation
• For each jet, require no other jet within \({\Delta}R\) = 1 (2.5 \(\times\) jet radial parameter) of the jet axis
• Methods exist for dealing with mislabelled training data
  • Majority vote filtering and consensus filtering can significantly improve classification accuracy for noise levels up to 30%\(^*\)

* Identifying Mislabeled Training Data, C. E. Brodley and M. A. Friedl, Journal of Artificial Intelligence Research, Vol. 11, pages 131-167, 1999

Data production

The data that I am using for the results shown are:

• \(pp\) at 7 TeV, high-\(p_\mathrm{T}\) data, Pythia 6, Perugia 2011 tune, full detector Monte Carlo simulation
• Charged-track jets reconstructed with Anti-\(k_\mathrm{T}\) with a radial parameter of 0.4
• Constituents: "Picotracks" with \(p_\mathrm{T}\) > 0.15 GeV
• Code exercised using the ALICE "EMCAL jet framework"
• No other cuts at this stage
• What about heavy-ions?
  • One step at a time...

Data preselection

• There were some jets for which a partonic flavour could not be assigned; these are discarded
• Heavy-flavour jets are discarded, as \(b\)-and \(c\)-tagging is not the topic of this analysis
• Single-track jets are removed. The justification for doing this is that many cases it is not possible to construct a meaningful observable for these jets
  • Probably not much we can do for these jets anyway, so they are removed
  • Try to recover these later maybe
• All remaining jets with NaN values were dropped; this amounts to less than 0.1% of jets (this might bias the data, but the impact will be extremely small)

Data munging

• Started with more than 300 features
• Removing duplicate and highly correlated features was critical for enabling my data to fit in RAM
  • I could reduce the size of my data by more than 50% with a conservative cut on correlations between features
  • Two-step process: cut on absolute value of Pearson correlation to identify linearly correlated features, then cut on absolute value of Spearman correlation to identify monotonically-related features
• Needed to reduce data size considerably (to fit into RAM) while retaining useful features
  • Succeeded in reducing number of features to ~50

Feature ranking & selection

• The question addressed by this work can be formally stated as follows: can we use quantitative characteristics of the jets to classify them as quark-jets or gluon-jets?
• This invites the question: how should we find the variables that provide the best discrimination between quark-jets and gluon-jets?
• We can use domain knowledge to drill down to what are believed to be the best discriminants; observables that:
  • Can explain most of the variance in the data
  • Are minimally correlated with each other
  • Provide the best predictive power
• How to optimally search the feature space? (Manual inspection may be impractical for a large feature set)

Fast feature selection

• The usual methods provided for feature selection work well with small-sized data, but become impracticably slow for data of the size that is typical in HEP
• It was necessary to devise a faster method
• I adapted a filter-based method that involves ranking by information gain (Kullback–Leibler divergence) or Gini impurity
• Basic idea: rank by one or both of these metrics, and do this repeatedly for a large number of bootstrap resamplings from the data and average
  • For each feature, we now have a nonparametric estimate of the mean and can calculate a confidence interval for that estimate!

Information gain & Gini impurity

• Information gain is based on the concept of entropy from information theory and is commonly used to decide which features to use when growing a decision tree

\[ Entropy = - \sum_{i}{p_i}{\log_{2}}{p_i} \]

• At each level, the feature with the highest information gain is chosen
• An alternative measure of "node impurity" commonly used in decision tree learning is the Gini impurity:

\[1 - \sum_{i}{p_i}^2\]

Information gain distributions

Note "BOGUS" random probes injected into the data

Feature ranking

Median information gain

Remove features with median information gain less than that for random probes, or within one standard deviation

Correlation analysis

• Filter-based methods for feature selection are fast, but they only look at each feature individually and do not take into account interactions between features
• Typically we would plot a correlation matrix to examine the Pearson (linear) correlations between features
• This visualisation type commonly used in particle physics, albeit with smaller matrices (\(<\) 20 features)
• "OK" for small matrices, but the information value of visualisation diminishes for larger matrices
  • Also: slightly more than half the plot is redundant!

Visualising a correlation matrix as a biplot

• Consider instead visualising correlation structures by performing principal component analysis (PCA) and plotting a biplot that shows pairwise correlations as angles:

\[ r\left(x_1,x_2\right) = \cos^2\left(x_1,x_2\right) \]

• In the following visualisation (next slide), PCA was performed on the data and the first two components extracted and plotted on a biplot (which also hints that existing features don't appear to be efficiently leveraging the apparent structure seen in the data)

Correlation matrix biplot

Screeplot

This screeplot shows that the most of the variance in the data, and hence most of the structure, is contained in perhaps no more than half a dozen principal components

Visualising a correlation matrix as a force-directed network

Connection strength between features (nodes) is proportional to their marginal correlation; the node colour provides a measure of feature importance

Visualising a correlation matrix as a force-directed network

One can also plot a partial correlation network that shows the conditional pairwise correlations between features

Visualising a correlation matrix as a force-directed network

We can remove edges between weakly-correlated predictors on a partial correlation network to isolate them

Prototype jet classifier

• I decided to use a boosted decision tree as my classifier
• Specifically, I used Stochastic Gradient Boosting or Gradient Boosting Machine (GBM)\(^\dagger\) in xgboost
• There are several justifications for this choice:
  • Speed: whereas R is single-threaded, XGBoost runs using multiple threads and is very fast\(^*\)
  • Outstanding performance in Kaggle competitions
  • Can do early stopping: quit training if performance degrades after N rounds; great for tuning\(^*\)
  • Can be tweaked to do Random Forests also\(^*\)
  • Robust with respect to class noise

\(\dagger\) Greedy Function Approximation: A Gradient Boosting Machine, Jerome Friedman, The Annals of Statistics, Vol. 29, No. 5, Oct., 2001; Additive logistic regression: a statistical view of boosting, Jerome Friedman, Trevor Hastie, Robert Tibshirani, The Annals of Statistics, 2000, Vol. 28, No. 2, 337-407

* I'm indebted to Szilárd Pafka for this info. He's done great work benchmarking these tools.

Results

Confusion Matrix and Statistics

          Reference
Prediction Gluon Quark
     Gluon 44991  8224
     Quark 34660 12125

               Accuracy : 0.5712          
                 95% CI : (0.5681, 0.5742)

• We correctly tagged 60% of true quark jets
• We mis-tagged 44% of true gluon jets ☹
• The main problem appears to be one of overfitting \(-\) performance on the training data is reasonable
• The task ahead is to deal with the overfitting problem to hammer down the mis-tag rate

Future work

• I have unbalanced classes: $>$80% of jets are gluon-jets
  • I got bad results until I realised that most of the classifiers I tried were optimising for accuracy; results improved when I optimised for AUC instead
• Down-sampling the majority class could hurt performance, as we remove points that could be useful for defining the optimal decision boundary
• Bagging should help ameliorate this
• I plan to CV bag\(^*\) to build an ensemble learner
• Add Random Forests and Deep Learning into the mix?

* Dissecting the Winning Solution of the HiggsML Challenge, Gábor Melis, Journal of Machine Learning Research, Workshops and Proceedings Vol 42, pp. 57–67, 2015

Summary

• It's often said that 80% of data analysis is spent on data munging\(^*\) \(-\) this was certainly true in my case
• However, I've found a good set of tools for streamlining this process
• I've shown how it is possible to use R for fast prototyping of a new method for binary classification
  • To the best of my knowledge, nobody has tried to do a particle physics analysis like this in R before
  • Will be invaluable for helping me decide where to direct my efforts later when building a final model

* Exploratory Data Mining and Data Cleaning, Dasu T, Johnson T (2003), John Wiley & Sons

Thanks!

alice-machine-learning@cern.ch