Apprentice for Event Generator Tuning

Apprentice is a tool developed for event generator tuning. It contains a range of conceptual improvements and extensions over the tuning tool Professor. Its core functionality remains the construction of a multivariate analytic surrogate model to computationally expensive Monte-Carlo event generator predictions. The surrogate model is used for numerical optimization in chi-square minimization and likelihood evaluation. Apprentice also introduces algorithms to automate the selection of observable weights to minimize the effect of mis-modeling in the event generators. We illustrate our improvements for the task of MC-generator tuning and limit setting.


Introduction
Monte Carlo-based (MC) event generators are necessary tools for interpreting data at the high energy frontier. MCs contain O(10-100) parameters that are tuned to match selected data to allow predictions for other sets of data. In many cases, the limiting factor to precision physics at the LHC is our lack of confidence in the MC predictions. The tuning task is daunting because of the large number of parameters that should be explored and the computational cost of simulations. For the large datasets of collider data that are available for MC tuning, a common heuristic to evaluate the goodness of a prediction is: where S O is the set of observables O used in the tune, each observable has a weight w O represented by a vector w, t b (p)/R b is the theory prediction/reference data in a given bin b of an observable and the ∆t b /∆R b 's are error estimates on these quantities. Our problem is to minimize χ 2 (p, w) as a function of the adjustable parameters p and possibly the observable weights w. To accomplish this, one needs a range of theory predictions t b for different possible parameter choices p and a method or principle for choosing w.
As an additional output, it is desirable to have an estimate of what range of parameters p are compatible with the data at a given confidence level.
In the following, we report on apprentice, a successor to the tool Professor [1] that was developed to accomplish these goals.

Problems Solved using apprentice
Professor is a public tool that was introduced to automate and accelerate the tuning of event generators. Based on a modest number of MC simulations, it constructs a polynomial approximation surrogate to these predictions in each bin of a histogram representing the observables that drive the tuning. Given the surrogate function, an "optimal" set of tuned parameters p is determined by minimizing the heuristic (1).
Professor has been used successfully to create a number of MC tunes and in other contexts. Nonetheless, there are certain aspects that could be improved. First, the polynomial fit to the MC predictions is adequate for well-behaved distributions, but does not perform well for a prediction x with distribution ∼ 1/x. Second, the suite of minimizers available is limited to those in SciPy and the CERN Minuit package. The dependence on minimizer, if any, needs to be explored. Third, the choice of weights w O is done manually, and any optimization over them must be done in a costly, iterative process. An algorithm for automatically selecting these weights is desirable.
The desire to address these three issues has led to the development of a new public tool that is sufficiently different from its predecessor to have its own name: apprentice.
In the following, we motivate the use of apprentice by listing the problems that apprentice can currently solve and by highlighting the advantages of using apprentice to solve these problems, especially for HEP applications.

Rational Approximation as a Surrogate Function
Polynomial models are relatively easy to build and use. However, they have poor extrapolation behavior and are severely limited in their ability to cope with singularities. These drawbacks can reduce their effectiveness at representing physics models. On the other hand, rational functions (quotients of polynomials) can be considerably more effective [2,3] at representing models that have real or apparent pole structure. Unfortunately, rational approximations can be numerically fragile to compute and are prone to having spurious singularities.
apprentice is capable of computing multivariate rational approximations r(x) = p(x)/q(x) to simulation data using three approaches. The first is based on the univariate methods of [4,5] and provides a robust and efficient way to compute the coefficients of p(x) and q(x). Although it tries to reduce the appearance of unwanted singularities by using ideas from linear algebra to minimize the degree of q(x), it does not guarantee that r(x) will be pole free in the parameter domain D.
The second uses a constrained optimization formulation that includes structural constraints on r(x) to enforce the absence of poles in D. Although this approach is computationally more expensive than the first, it guarantees that the computed approximation is free of poles for box-shaped parameter domains, which can be crucial when computing surrogate models for use in optimization. In particular, the guaranteed absence of poles ensures that subsequent optimization problems involving our rational approximations are well-defined.
The third approach is a specialization of the second that is used to efficiently find a polefree rational approximation by making q(x) = A T x+b linear. This is achieved using auxiliary variables to write the dual of the linear program min x A T x + b s.t. x ∈ [L, U] as a constraint to force the denominator to be greater than 0 or be pole-free. The advantage of this approach is that it requires solving the optimization problem only to find a pole-free rational approximation (as opposed to an iterative approach) and hence it is more efficient.

Tuning Problem
The goal of the tuning problem is to find an optimal set of physics parameters, p * , that minimizes the function (1). When the theory prediction t b (p) is based on simulated data from the MC event generator, the computational cost is high (the generation of 1 million events for a given set of parameters consumes about 800 CPU minutes on a typical computing cluster), severely limiting the number of parameter choices p that can be explored in the tuning. To overcome this issue, we construct a parametrization of a modest number of MC simulations i.e., of MC b (p) and ∆MC b (p) of each bin b as a polynomial or rational approximation using apprentice as discussed in the previous section. Thus t b (p) becomes our analytic surrogate f b (p), which is easy to compute (it can be evaluated in milliseconds) and minimize. For performing this optimization, the solver options of Truncated Newton method [6], Newtonconjugate gradient method [7], LBFGS-B algorithm [8], and Trust region algorithm [9] from python's SciPy package are provided within apprentice.

Automatic Selection of Observable Weights
To tune parameters as discussed in the previous section, the weights w O in χ 2 (p, w) need to be specified. In practice, the weights are adjusted manually, based on experience and physics intuition: the expert fixes the weights and minimizes (1) over the parameters p. If the resulting fit is unsatisfactory, a new set of weights is selected, and the optimization over p is repeated until the tuner is satisfied. Our goal is to automate the weight adjustment, yielding a less subjective and less time-consuming process.
In apprentice, the automatic selection of weights is achieved using two mathematical formulations: bilevel and robust optimization. In bilevel optimization, a merit function is minimized over the set of weights. apprentice provides three options for merit functions. The first is the portfolio objective function, which is motivated by portfolio optimization in finance, where the goal is to maximize the expected return while minimizing the risk. Translated to our problem, we want to minimize the expected error over all observables while also minimizing the variance over these errors. The second and the third merit functions are based on the central values for scoring schemes with a scoring rule, S (P, where model P has mean performance µ P and variance σ 2 P . For our application, x corresponds to the simulation prediction f b (p), µ P to our observation data R b , and the variance σ 2 P to our data uncertainty ∆R b . The second merit function maximizes the sum over all observables of the mean across all bins in that observable of each bin's scoring function. The third merit function maximizes over the median. The scoring function for bin b is defined Since we don't have the full analytical expression of the outer objective function, this function is approximated with a radial basis function (RBF) [10] to perform the bilevel optimization in apprentice.
Robust optimization is a single-level formulation for finding the optimal weights for χ 2 (p, w). It estimates the parameters p that minimize the largest deviation ( f b (p) − R b ) 2 over all bins in an uncertainty set U b of bin b. Assuming that the experiment and the MC simulation are described using independent random variables with mean R b , the uncertainty set U b for each bin b is described by the interval 2 is a square function, the largest deviation occurs at either end of the interval. Using this, the robust optimization formulation can be easily rewritten as a one shot single level optimization problem. To avoid the trivial solution of all weights being zero, the sum of the weights w ∈ [0, 1] across all the observables is required to cover at least µ% of the total observables in S O , where µ is a hyperparameter set by the end user.

Using apprentice
The source code for apprentice is available at https://github.com/HEPonHPC/apprentice. For each problem described in Section 2, we provide a convenient script that the end user can use directly with appropriate input arguments to solve that problem for their use case. In this section, we first describe how to clone and install apprentice, and, then, how the interface scripts can be used by the end user to directly solve their problem. For this, we describe the important arguments of each script in this section. For the complete usage, run the script using the -h flag. Then, we describe what output to expect from each of these scripts. Additionally, we also note the software features that are currently implemented for each problem.

Creating Polynomial and Rational Approximations
We provide an easy-to-use script at apprentice/bin/app-build to create a rational approximation r(x) = p(x)/q(x) for an arbitrary degree for p(x) and q(x) or when the order of q(x) is one or zero i.e., polynomial approximation. The important arguments of the script are:

Solving the χ 2 Tuning Problem
For solving the tuning problem described in Section 2.2, the interface script is located at apprentice/bin/app-tune2. The output of this script are the parameters p * , χ 2 (p * , w), and some other pertinent information about the optimization for logging purposes. This script also supports running the optimization with multiple starting points in parallel using mpi4py [11].

Event Generator Tuning by Automatic Selection of Observable Weights
In this section, we describe the interface scripts for solving generator tuning problems with automatic selection of observable weights.

Bilevel Optimization Formulation
Solving the generator tuning problem with selection of observable weights using the bilevel optimization formulation as described in Section 2.3 requires two steps: (a) generating the data to train the RBF that will be used as an approximation of the outer level objective function and (b) solving the bilevel optimization problem with either of the three merit functions of the portfolio function, the mean of the scoring function (meanscore), or the median of the scoring function (medianscore). The interface script for generating the training data is at apprentice/pyoo/bin/pyoo-train. The arguments args[0], args [1], args [2], -e, -s|survey, and -r|-restart of this script mean the same as those described in Section 3.3. Additional (important) arguments of the script are: (a) -d|-design Size of the initial design, (b) -o|-output Training data output file, and (c) -seed Random seed to use.
Using the training output data (a JSON file), we can run the interface scripts for solving the bilevel optimization using the portfolio objective found at apprentice/pyoo/bin/pyoo-run-portfolio, meanscore objective found at apprentice/pyoo/bin/pyoo-run-meanscore, and medianscore objective found at apprentice/pyoo/bin/pyoo-run-medianscore. The arguments for all three scripts are the same. Also, the arguments args[0], args [1], args [2], -e, and -s|-survey of these scripts mean the same as those described in Section 3.3. Additional (important) arguments of the script are: (a) args [3]: Training data JSON file, (b) -n|-niterations: Number of iterations to use, and (c) -o|-output: Output JSON file.
The output of each script is a JSON object with keys such as X_outer, X_inner, Y_outer, hnames and pnames. The keys X_outer and X_inner point to an array of array where each array element in the outer array are the weights w * of the observables and the parameters p * obtained in each iteration of the outer optimization. The best weight and parameter element in the output corresponds to the index of the smallest value in the array against the key Y_outer i.e., smallest outer objective value. The order of the observables and parameter dimensions that the weights and parameter values correspond to can be found against the keys hnames and pnames, respectively.

Single Level Robust Optimization Formulation
To solve the generator tuning problem where the observable weights are selected using the robust optimization approach (see Section 2.3), the interface script is located at apprentice/pyoo/bin/pyoo-robopt. The important arguments of the script are: (a) -w|weights: Weight file (see description of arg[0] in Section 3.3), (b) -d|-expdata: Measurement data (see description of arg [1] in Section 3.3), (c) -a|-approx: MC central values approximation JSON file (see description of arg [2] in Section 3.3), (d) -e|-error: See description of -e|-error in Section 3.3, (e) -o|-outtdir: Output directory where the output JSON file and the optimization log files are stored (more details below), (f) -s|-solver: Solver where the options are ipopt that makes use to IPOPT solver with AMPL solver interface (ASL) using PYOMO [12,13] and cyipopt [14] uses the python wrapper around IPOPT and avoids making use of PYOMO and the ASL, and (g) -m|-mu: Hyperparameter µ ∈ [0, 1] is the value to be used in the constraint to avoid tivial solutions of all weights being zero.
The output of the script is a JSON object where the weights w * of the observables are in the array pointed by the key X_outer, the parameters p * , which is obtained by solving the χ 2 tuning problem using w * , are in the array pointed by the key X_inner. The order of the observables and parameter dimensions that the weights and parameter values correspond to can be found against the keys hnames and pnames, respectively. Additionally, the object pointed by the key log contains some pertinent information about the optimization for logging purposes.

Impact of apprentice for HEP applications
In this section, we discuss the impact of the problems that can be solved using apprentice for HEP applications.

Rational Approximation for High Energy Physics
As a demonstration of the utility of rational approximations, we consider the problem of setting limits on dark matter particles in a future direct detection Xenon-based experiment (see [15,Sec. 26]). The physics model has three parameters, x = (m χ , c + , c π ), representing the dark matter particle mass (ranging over 10-100GeV/c 2 ) and two (dimensionless) couplings to ordinary matter (ranging from 10 −4 -10 −3 and 10 −3 -10 −1 , respectively). We consider a "binned likelihood" analysis [16] of L(x (k) |d 1 , d 2 , . . . , , where the d b represent six data points and N b (x (k) ) denote the simulated quantities for a point x (k) that correspond to the d b . Since no measured data exists as yet, we assume a signal consistent with a dark matter mass m χ = 10 GeV/c 2 and an interaction strength large enough to produce approximately 100 events, yielding {d 1 , d 2 , . . . , d 6 } = {70.4, 26.7, 9.8, 3.4, 1.0, 0.2}, simulated by using a specific framework of the generalized spin-independent response to dark matter in direct detection experiments [17,18]. Regions of parameter space consistent with the simulated data are found using MultiNest [19][20][21]. This requires the evaluation of the likelihood function at tens of thousands 2 of x (k) to succeed, and the computational cost can be substantial. To reduce the cost, we replace the expensive simulation of N b with the cheaper evaluation of a rational approximation surrogate. The results presented here are for rational approximations of degree M = 4, N = 4 as well as polynomial approximations of degree 7. 3 The maximization of the likelihood function requires about 30,000 evaluations in all cases, but is about a factor of 50 times faster using the rational or polynomial approximation. Figure (1) shows two-dimensional profile-likelihood projections of the three-dimensional likelihood limits when using the expensive full simulation, our rational approximation results, and a polynomial approximation. Clearly, the rational approximation faithfully represents the full simulation, while the polynomial approximation is only valid over a limited range.  Figure 1: Two-dimensional profile-likelihood projections of a 3-dimensional parameter space with superplot [22]. Regions of higher likelihood are darker. The data are normalized to the maximum observed likelihood.

Tuning HEP Event Generators by Automatic Selection of Observable Weights
As described earlier, we have developed several algorithms for adjusting the weights used in our tuning heuristic. For comparison, we re-consider a tuning exercise performed by the ATLAS collaboration [23, 24] using a large set of LHC data. Figure 2 shows cumulative distribution of χ 2 values per bin for several different classes of observables. The results of the ATLAS expert tune and a much simpler approach using equal weights for all observables is also shown. Note that the parameters obtained from the robust optimization perform well for Jet shapes and Track-jet UE. We also see that near the variance boundary, the parameters obtained from the Expert tune perform better for Multijets and tt gap whereas the parameters obtained from the other approaches perform better for Substructure. These plots demonstrate the trade-off in fitting among the different approaches, which enables the physicist to use these results as guidance for selecting the most appropriate tuning method depending on the categories that are of greater significance.

Conclusion
In this paper, we propose an improved software package called apprentice that enables the construction of pole-free rational approximation, allows for solving the χ 2 minimization problem using a various optimization algorithms that can detect multiple local minima, and provides multiple optimization formulations to automatically assign weights to the observables yielding a less subjective and less time-consuming process to find the optimal event generator tune. In the future, we are interested in (a) tuning event generators using the MC generator data directly using derivative free optimization methods, (b) providing non-linear optimization algorithms to minimize χ 2 that are robust against all local minima, and (c) exploring machine learning techniques for dimensionality reduction and for improving the automatic weight assignment toward obtaining optimal generator tune. If your are interested in these problems or in any other interesting extensions to apprentice, we look forward to hearing from you. To contribute, either raise an issue in the apprentice project or contribute via code by forking the project. The apprentice project is located at https://github.com/HEPonHPC/apprentice.