Reducing the dependence of the neural network function to systematic uncertainties in the input space
Stefan Wunsch Simon Jörger Roger Wolf Günter Quast
physics.data-an












Stefan Wunsch1,2 (corresponding author) stefan.wunsch@cern.ch
Simon Jörger1 simon.joerger@cern.ch
Roger Wolf1 roger.wolf@cern.ch
Günter Quast1 guenter.quast@kit.edu
1 Karlsruhe Institute of Technology, Institute of Experimental Particle Physics, Karlsruhe, Germany

2 CERN, Geneva, Switzerland




Reducing the dependence of the neural network function to systematic uncertainties in the input space
    Stefan Wunsch Simon Jörger Roger Wolf Günter Quast
    
  First Draft: July 26, 2019 

  Current Draft: August 31, 2020
=====================================================================================================





Applications of neural networks to data analyses in natural sciences are complicated by the fact that many inputs are subject to systematic uncertainties. To control the dependence of the neural network function to variations of the input space within these systematic uncertainties, several methods have been proposed. In this work, we propose a new approach of training the neural network by introducing penalties on the variation of the neural network output directly in the loss function. This is achieved at the cost of only a small number of additional hyperparameters. It can also be pursued by treating all systematic variations in the form of sta­tis­ti­cal weights. The proposed method is demonstrated with a simple example, based on pseudo-ex­pe­ri­ments, and by a more complex example from high-energy particle physics.








§ INTRODUCTION



Neural network (NN) techniques are in wide and increasing use to solve classification and regression tasks in the analysis of high-energy particle physics data. Examples of their use in physics object identification, e.g. at the LHC experiments ATLAS and CMS, are the classification of particle jets induced by heavy flavor quarks <cit.> and the identification of τ leptons <cit.>. Examples for data analyses that make use of NNs not only for object identification, but to distinguish between signal- and background-like samples are the latest analyses of Higgs boson events in association with third generation fermions, at the LHC <cit.>. These classification tasks usually aim at the distinction of a signal from one or more background processes. They are characterized by a relatively small number of input parameters to the NN, of one or two orders of magnitude, which may reveal non-trivial correlations among each other.

Each physics measurement is subject to systematic uncertainties, which have to be propagated from the input space 𝐱={x_i} to the NN output f(𝐱). This usually happens in terms of variations of a given input parameter x_i within its uncertainties Δ_i. We abbreviate the set of Δ_i by Δ={Δ_i} and the set of modified input parameters by 𝐱+Δ={x_i+Δ_i}. These variations may be implemented in the form of variations of the actual values of x_i, or such that a sample, with a given value of x_i, enters the analysis with a different statistical weight, also referred to as reweighting throughout this text. Unlike varying the values of x_i, reweighting does not rely on a reprocessing of the dataset and therefore generally implies significantly smaller computational costs.

The possibility to implement prior information about systematic uncertainties already in the NN training is motivated by two considerations: Firstly, a powerful distinction between classes in principle, can be considerably compromised by systematic uncertainties. Integrating prior knowledge of uncertainties in the NN training helps in guiding the NN to focus on features in the input space that are less prone to such a performance degradation. This may even result in a gain for the analysis performance, as observed in Ref. <cit.>.
Secondly, the dependence of a systematic variation of a given feature x_i on other parameters {x_j}, j≠ i in the input space, might only be poorly known, or even unknown, and the user might want to generally uncorrelate the NN output from this uncertainty to assure a reliable response of the NN to the given task. Both points raise interest in training the NN with the boundary condition that the dependence of f(𝐱+Δ) on Δ should be minimal.

One way of achieving this decorrelation of f(𝐱+Δ) from Δ that has been proposed in the past and that we will refer to in more detail throughout this paper, makes use of a secondary NN that is trained in addition to the primary NN in an iterative procedure, resulting in an adversarial architecture <cit.> for robust binary classification <cit.>. This secondary NN has the task of drawing information of the systematic variation from the output of the primary NN. The output of the secondary NN is then included in the loss function of the primary NN as part of a minimax optimization problem. The resulting setup becomes insensitive to the systematic variation of the inputs. This method requires a relatively complex iterative training procedure; it introduces a large and to some extent arbitrary number of new hyperparameters implied by the choice of the architecture of the secondary NN, and requires the resampling of x_i within its uncertainties Δ_i.

Another approach to decorrelate f(𝐱+Δ) from Δ is to include the knowledge about systematic uncertainties in a systematics-aware objective function as proposed in Refs. <cit.> and <cit.>. An approach related to boosted decision trees is implemented by splitting the tree nodes using the signal significance including systematic uncertainties as objective, resulting in a classifier that successfully reduces the impact of systematic uncertainties on the result <cit.>. A similar approach for NNs has been studied in Ref. <cit.>. A comparison of systematics-aware learning techniques in high-energy particle physics has been carried out in Ref. <cit.>. In addition to the adversarial approach discussed above this study includes a comparison to data perturbation and augmentation, and tangent propagation <cit.>.

In our approach we implement a penalty on the differences between the NN output obtained from the nominal value of x_i and its variations Δ_i, directly into the loss function. For this purpose we use histograms of f(𝐱) and f(𝐱+Δ) filled during each training batch. The number n_k of histogram bins {k}, and the batch size n_b are hyperparameters of the training. To guarantee a differentiable loss function for the optimization of the trainable parameters of the NN, the histogram bins are blurred by a filter function applied to each sample b of the training batch, affected by the uncertainty variations, where b corresponds to a single sample represented by a point associated to each respective training dataset in the input space 𝐱. We use Gaussian functions 𝒢_k(𝐱), normalized to max(𝒢_k(𝐱))=1 as filters, where the mean and standard deviation are given by the center and half-width of histogram bin k. The count estimate can then be written as 𝒩_k(f(𝐱)) = ∑_b𝒢_k(f(x_b)), and the loss function consists of the two parts

    L_Λ =  L^' + λ Λ(𝐱, Δ) 
        
    
    
    
           with:
       Λ(𝐱,Δ) = 1/n_k∑_k(𝒩_k(f(𝐱)) - 𝒩_k(f(𝐱+Δ))/𝒩_k(f(𝐱)))^2,

where L^' corresponds to the loss function of the primary task, like for example the cross-entropy function for a classification task, and Λ(𝐱, Δ) to the term that penalizes differences in the NN function between f(𝐱) and f(𝐱+Δ). The factor λ controls the influence of the penalty and adds another hyperparameter to the training. The count estimate 𝒩_k(f(𝐱+Δ)) can be derived from 𝒩_k(f(𝐱)) in terms of reweighting, such that no reprocessing of the dataset during the training procedure is required.

In this approach more than one uncorrelated uncertainty simply adds to the sum of Λ_ℓ(𝐱, Δ), for ℓ uncorrelated uncertainties. Two fully (anti-) correlated uncertainties should be represented by a common variation for both uncertainties at the same time. While an exact modeling of correlations across uncertainties may not always be exactly known this knowledge is not strictly required by the method, as long as the loss function converges to its minimum and solves the defined task. The parameters λ_ℓ correspond to further hyperparameters, whose values relative to each other define different tasks of the NN training. We would like to emphasize that the use of a histogram of f(𝐱) (and f(𝐱+Δ) respectively) in the loss function might lead to a suboptimal performance with respect to the direct use of f(𝐱). Also we do not claim the resulting discriminator to be optimal for the final measurement.

In Section <ref> we demonstrate the method on a simple example based on pseudo-experiments. A more complex analysis task typical for high-energy particle physics is studied in Section <ref>. We summarize our findings in Section <ref>.



§ APPLICATION TO A SIMPLE EXAMPLE BASED ON PSEUDO-EXPERIMENTS



To illustrate our approach, we refer to a simple example based on pseudo-experiments that has also been used in Ref. <cit.>. It consists of two variables x_1 and x_2, which are the input to separate two classes, in the following labeled as signal and background. The input space is visualized in Fig. <ref>. A systematic uncertainty for the background class is introduced by two variations of x_2 by ±1. We consider only the discrete variations that quantify the difference between 𝒩_k(f(𝐱)) and 𝒩_k(f(𝐱+Δ)), which is sufficient to define the part of Λ(𝐱,Δ) to be minimized during the training process. We have checked that a Gaussian sampling with a standard deviation of σ=1, as applied in <cit.> would lead to the same result in a more complex setup.


([ 0 0 ])

([ 1 1 ])

([    1 -0.5; -0.5    1 ])

([ 1 0; 0 1 ])


The NN used to solve the classification tasks consists of two hidden layers with 200 nodes each, with rectified linear units as activation functions <cit.> and a sigmoid activation function for the output layer. The trainable parameters are initialized using the Glorot algorithm <cit.>. The optimization is performed using the Adam algorithm <cit.> with a batch size of e3. Our choice for L^' is the cross-entropy function. For Λ, we use 10 equidistant bins in the range [0, 1] of the NN output. We have not observed any significant performance dif­fe­ren­ces by varying the number of histogram bins within reasonable boundaries, though. Finally we set λ to 20. The training on 5e4 events is stopped if the loss obtained from the training dataset has not decreased for five epochs in sequence, on an independent validation dataset of the same size. In addition, we use e5 events for testing and to produce the figures to illustrate the result. The impact of the systematic variations on the NN output is shown in Fig. <ref> for the case of a classifier trained with a loss function given only by L^' (f_L^') and a classifier based on a loss function including the additional penalty term Λ (f_L_Λ).



As can be seen from Fig. <ref>, the approach successfully mitigates the dependence of the NN output on the variation of x_2 and therefore results in a classifier that is more robust in the presence of this systematic uncertainty. This is achieved on the expense of obliterating at least parts, if not all, separating information of x_2. Fig. <ref> visualizes the NN output as a function of the input space spanned by x_1 and x_2. The additional penalty term, Λ, leads to the intended alignment of the surface of the NN output with the variation of x_2, resulting in similar values of the NN output for all realisations of the systematic variation. We find our approach to have an effect similar to the setup described in <cit.>.





§ APPLICATION TO A MORE COMPLEX ANALYSIS TASK TYPICAL FOR HIGH-ENERGY PARTICLE PHYSICS



In the following, we apply the proposed method to a more complex task typical for high-energy par­ti­cle phy­sics. We use a dataset that has been released for the Higgs boson machine learning challenge described in Ref. <cit.>. This challenge uses a simplified synthetic da­ta­set from simulated collisions of high-energy proton beams with underlying hypothesized signal and background processes at the CERN LHC. The original target of the challenge was to separate events containing the decay of a Higgs boson into two tau leptons (signal) from all other events (background), to serve as benchmark for the success of different machine learning algorithms. The consideration of uncertainties, as required for a complete analysis of the data was not part of it. The dataset contains 30 input parameters, whose exact physical meanings are given in Ref. <cit.>. We split the dataset and use one third for training and validation of the NN and two thirds for deriving the following results.

For our example, we use all parameters as input for the NN training. In addition, we introduce a systematic uncertainty, resembling the fact that the momentum and energy of a particle are the results of external measurements with a finite resolution. For our study we assume an uncertainty of ±3% <cit.> on the transverse momentum of the reconstructed hadronic τ decay p_t^τ, measured in GeV and labeled as  in Ref. <cit.>. The distributions of the nominal and varied input parameters are visualized in Fig. <ref> (upper row). To allow for migrations in and out of the selected input space due to the systematic variation we restrict the originally available dataset by raising the lower p_t^τ requirement from 20 to 25 GeV. For the background the distribution of p_t^τ is steeply falling. Thus the variation is dominated by migration effects at the lower p_t^τ boundary, resulting in an overall normalization uncertainty. The signal shows a maximum around p_t^τ≈25 GeV, leading to a more apparent additional variation of the shape of the p_t^τ distribution, as shown in Fig. <ref> upper right. The dataset used for these results contains 814.9 (163750) weighted (unweighted) signal events and 162705.0 (238778) weighted (unweighted) background events using an additional scaling of the weighted number of signal events by a factor of two.

Instead of resampling the signal and background datasets with the varied values of p_t^τ, we introduce the systematic variation in form of statistical weights. In this way we give a higher (lower) statistical weight to subsamples with low (high) values of p_t^τ with respect to the nominal sample. These weights are determined from the p_t^τ distributions shown in Fig. <ref> (upper row) for the background and signal sample, respectively. By construction all correlations across features of the input space are conserved by the reweighting, thus that reweighting p_t^τ leads to shape variations also of correlated observables, e.g., like the reconstructed missing transverse momentum or the estimate of the invariant di-τ mass, described in Ref. <cit.>, as shown in Fig. <ref> (lower row). We would like to emphasize that this reweighting technique is in fact the only way to apply a systematic variation of p_t^τ that respects the correlations to all other features of the input space on the given dataset. In a realistic analysis the reweighting technique is not meant to replace the resampling, but rather to complement it. A resampling could and should be applied, where correlations across input features may not be desired. To give an example, p_t^τ is mostly determined from track information. Therefore an uncertainty in the missing transverse momentum due to uniformity uncertainties in the calibration of the hadronic calorimeter should not impact p_t^τ with a correlation of 100%. As in the case of the simple example of Section <ref> we use only the two discrete shapes corresponding to the ±3% shifts in p_t^τ, which are a sufficient input for the minimization of the loss function during the training process. Samples of intermediate realizations of these shifts have been checked to lead to the same result despite the more complex setup.



The NN has the same architecture as described in Section <ref>. For the implementation of f_L_Λ we chose 20 equidistant bins in the range of [0, 1] of the NN output for Λ, and λ=20. The batch size is set to 10^3. The optimization of the trainable parameters is performed on 75 of the training dataset and stopped if the loss has not decreased for 10 epochs in sequence, on the remaining part of the training dataset. The results are shown on an independent test dataset. We would like to emphasize that λ=20, is a free choice that has been made for illustrative purposes only. In a realistic application the optimal choice of λ should be studied on a case by case basis.



In Fig. <ref> the NN outputs f_L^' and f_L_Λ are shown. As in the case of the simple example given in Section <ref>, though less pronounced, the training based on a loss function including Λ leads to a mitigated dependence of the NN output on the systematic variation of p_t^τ. An important difference between both examples is that the uncertainty of the simple example given in Section <ref> is exclusively shape altering. In contrast to this the uncertainty variation in this more complex example includes a significant component acting on the normalization of the NN output, especially for the background distribution. A pure normalization uncertainty that does not lead to noticeable differences in the input space that can be related to its systematic variation can not be mitigated. In consequence a dominant overall normalization uncertainty, visible especially for the background distribution of f_L^', is not significantly reduced by the use of f_L_Λ.

In Fig. <ref> the p_t^τ distributions for signal and background for the full sample, and for two signal-enriched subsamples are shown. The latter are obtained by a restriction of f_L^' and f_L_Λ to a value larger than 0.7. On the full sample a generally harder p_t^τ spectrum for the signal is observed with a maximum around 45 GeV, in contrast to a steadily falling and softer spectrum for the background. In the signal-enriched subsample based on f_L^'>0.7 the p_t^τ distribution for the background is biased towards the same distribution as for signal. In the signal-enriched subsample based on f_L_Λ>0.7 this bias is alleviated and the p_t^τ distributions for signal and background are qualitatively unchanged with respect to the full sample.



At the LHC experiments the presence of the Higgs boson signal has been inferred from hypothesis tests based on a likelihood ratio between the case of including the Higgs boson signal and that of the null hypothesis without Higgs boson signal <cit.>. Systematic uncertainties have been incorporated in form of nuisance parameters, which might be correlated, e.g., across processes, into the likelihoods. Best estimates and constraints on these nuisance parameters have been obtained by nuisance parameter optimization. The presence of the signal has been quantified, e.g., by means of its statistical significance in terms of Gaussian standard deviations (s.d.), in the limit of large numbers. To serve our discussion we emulate this discovery scenario, in a simplified way, constructing binned likelihoods for the signal and null hypotheses based on the histograms shown in Fig. <ref>. In addition to the statistical uncertainties of the pseudo-data we incorporate the uncertainty indicated by the bands in Fig. <ref> as process- and bin-correlated variations in the likelihoods, bound to a single nuisance parameter θ, following the prescriptions of <cit.>. The fit of a Higgs boson signal hypothesis with a single signal strength parameter of interest, μ, to the pseudo-data, including the signal as expected by theory, leads to a constraint of the uncertainty in θ to 3% of its initial value, both in the case of f_L^' and f_L_Λ as input distributions to the fit. This constraint is dominated by the power of the pseudo-data to determine the normalization related to θ, especially in the first bins of the background dominated pseudo-data sample distribution, e.g., with more than 65 thousand counts in the first bin. When splitting the uncertainty into two independent nuisance parameters, θ_norm to govern the pure normalization uncertainty, and θ_shape to govern the pure shape altering uncertainty, we find the initial normalization uncertainty to be 7.6% (2.2%) for the background (signal) sample. We anticipate that the implementation with two independent nuisance parameters is not fully correct, but keeping this caveat in mind the study still serves the test we are interested in. After the fit of the Higgs signal hypothesis to the pseudo-data we observe the same constraint as on the uncertainty in θ before on the uncertainty in θ_norm. We observe an ≈35% correlation between θ_norm and μ. The constraint on the uncertainty in θ_shape is 0.8 (0.4) for f_L^' and f_L_Λ as input distributions to the fit, with a correlation of 55% (5%) to μ. We observe similar results when performing a fit of the null hypothesis. The reduction of the correlation of θ_shape with μ, when using f_L_Λ instead of f_L^' gives a quantitative measure in this case of the decorrelation of the shape altering part of the uncertainty with the parameter of interest.

In Fig. <ref> the significance of the analyzed signal in the pseudo-data, based on the fit to the null hypothesis is shown as a function of the hyperparameter λ, where λ=0 corresponds to f_L^' as input to the fit. Using f_L^' as input to the fit leads to a significance of 6.7 s.d., corresponding to a combined systematic and statistical relative uncertainty in the parameter of interest of Δμ/μ=15%. This significance drops to a value of 5.2 s.d., corresponding to Δμ/μ=19%, for λ=20. Such a drop is expected, since p_t^τ plays an important role in the separation of signal and background, not only as a single feature, but also via its correlations to other features in the input space <cit.>. The scan of λ in this way visualizes to what extend the separation relevant information related to p_t^τ in the input space that is vulnerable to the variation of p_t^τ, is masked during the training process for increasing values of λ. The information loss seems small for values of λ≤ 5 with a significant drop around λ≈10 and a plateau around λ≈20, which is the value we have chosen for our study. At this point most of the separation relevant information related to p_t^τ that is vulnerable to the variation of p_t^τ seems to be masked out from the training, such that f_L_Λ turns mostly blind for p_t^τ. Implicitly this can also be inferred from Fig. <ref>, where the distribution of p_t^τ qualitatively is the same for the signal-enriched and the inclusive samples.

In turn the uncertainty on the significance due to the systematic variation drops, roughly proportional to the loss in significance, from ≈7.5% (for λ=0) to ≈1.8% (for λ=20). We estimate the contribution of the systematic variation in p_t^τ to Δμ/μ, with 6.6% (for λ=0), dropping to 1.8% (for λ=20). At the same time, and with a larger slope, the absolute contribution of the statistical uncertainty to Δμ/μ increases from 13.4% (for λ=0) to 19.3% (for λ=20), resulting in the overall decrease of the significance for increasing values of λ, for the given example. The loss in statistical power stems from the worse separation of signal and background for increasing values of λ, as also visible from Fig. <ref>.

Increasing λ to larger and larger values leads to another drop of the significance, which converges to the value for a single counting experiment that does not distinguish between signal and background, in the limit of λ→∞. This can be understood in terms of Λ completely dominating the loss function thus that L^' will more and more loose influence in the training task. As a consequence the NN will primarily be optimized on the suppression of the variation of p_t^τ rather than the separation of signal and background.

We would like to point out at the end of this discussion that it is usual practice in a measurement scenario to accept the increase of statistical uncertainty, which can in principle be controlled by an increase of the dataset for the benefit of a reduced sensitivity of the measurement on systematic variations of its input parameters, which might be difficult to control. We anticipate though that in the given scenario λ=0 remains the choice that maximizes the significance of the analysis despite its larger sensitivity to the systematic variation in this case. Our choice of λ=20 should be viewed as a free while still sensible choice to showcase the reduction of the influence of the systematic variation on the NN output.





§ SUMMARY



We have presented a new approach to reduce the dependence of the NN output to variations of features x_i of the NN input space due to systematic uncertainties in the measured input parameters. We achieve this reduction by including the variation of the NN output w.r.t. the nominal value of x_i in the loss function used for training. Compared to a previously published method of using an adversarial technique, the complexity of the presented method is reduced to one additional term in the loss function with less hyperparameters and no further trainable parameters. Systematic variations can be inscribed in the form of statistical weights, implying no further needs of reprocessing, further reducing the complexity of the training. Additional uncertainties just add to the sum of penalty terms in the loss function. In turn the method requires batch sizes large enough to populate the blurred histogram of the NN output used for the evaluation of the variation w.r.t the nominal value of x_i in the loss function.

We have demonstrated the new approach with a simple example directly comparable to a solution of the same task exploiting the adversarial technique, and a more complex analysis task typical for high-energy particle physics experiments. In all cases the dependence of the NN output on the variation of a chosen input parameter is successfully mitigated. In application to a high-energy particle physics measurement this leads to a result less prone to systematic uncertainties, which is of increasing interest in the presence of growing datasets, where statistical uncertainties play a subdominant role in the measurement.





splncs


