Dimensionally Constrained Symbolic Regression

1 Introduction

Extraction of a physical parameter from data is an area where application of symbolic regression can be useful, especially if the relationship between the parameter of interest and measured variables are contrived and non-linear [2]. In the problem of mass measurement of a particle in high-energy physics, the relationship is determined exactly if all the decayed particles are detected and measured by the detector. The problem becomes non-trivial if some of the particles escape detection.

In $W$ boson mass measurement with $W \to ℓ ν$ at hadron colliders, the neutrino ( $ν$ ) escapes detection, but its transverse components of the momentum are measured indirectly. In this case, transverse mass ( $M_{T}$ ) is known to be the most sensitive variable to the $W$ boson mass $M_{W}$ [1].

In searches of interest at LHC and the Tevatron, such as the Higgs particle and supersymmetry phenomena, there are two or more particles which escape the detector without leaving any signal. It makes mass measurement of these particles challenging. In some of cases, there many solutions are known, but it is not clear in what sense they are optimal for the mass measurement [4, 5, 3]. By applying symbolic regression to these problems, we may be able to derive optimal equations sensitive to mass and gain new insight.

2 Dimensionally Constrained Symbolic Regression and Its Implementation

Symbolic regression(SR) is a method which can be used to derive relationship between input variables and the desired output values. It is an application of genetic programming technique. In SR, an individual of a population is symbolic form of an equation. By applying genetic operations, the population evolves and gradual improvement of the overall population as well its most fit individual occur. Selection is done at the individual level at each generation based on its fitness.

We need to define the representation, functions/operators to be used, terminals that will be used to form an individual in SR. These, in addition to genetic operations, fitness function, and numerous parameters, define an SR. Later, we will describe the dimensionally constrained symbolic regression (DCSR) used for mass reconstruction problems in HEP. To have the flexibility DCSR as an option, we implemented a general SR in C++ language based on ROOT C++ libraries [7].

2.1 Representation and Genetic Operations in SR

Internally, mathematical expression is represented as a binary tree (Fig. 1). Operators or functions appear at tree nodes and arguments to them appear as branches. Binary tree representation allows for efficient manipulation and evaluation of mathematical expressions.

The functions and operators allowed are as follows:

Basic arithmetic operators: $+$ , $-$ , $\times$ , $\div$
Negation
Power
Transcendental functions: $\sqrt{}$ , $log$ , $sin$ , $cos$ , $exp$
Unit step function

For a function or operator that accepts only one argument, the second leaf may be populated, but is not evaluated. Although the second argument has no meaning with respect to the function or operator, it may nonetheless affect evolution. With cross-over operations and mutations, this second leaf may be picked up for genetic operations. Our implementation is flexible enough such that we can add more functions or operators, if needed.

Figure 1: A binary tree of expression $(v 1 + v 2) * v 3$ .

When creating an arbitrary expression, the functions or operators are picked randomly according to certain probabilities. The probabilities should be provided as an input to SR. To a certain degree, they will depend on the problem to be solved. We can control the relative rates these functions or operators are chosen through the input parameters of SR.

2.2 Terminals

The input variables and constant become terminals of expression trees. In our implementation of SR, the relative rate for variables and constants can be controlled. However, we set the probability to choose among the variables the same. For the DCSR, which we will describe latter section, we only allow dimensionless arbitrary constants.

2.3 Building Expressions from Functions/Operators and Terminals

A random expression can be build by first choosing a function or an operator randomly with the predefined probability. This becomes the head node of a tree. Then the leaf nodes are built recursively by creating a new expression. An expression can be any of the following:

Terminal - constant or variable
Unary function - $f (x)$ , where $x$ is another expression
Binary function - $g (x, y)$ , where $x$ and $y$ are expressions

Since the expression can go on forever, probability to choose a terminal must not be 0, and/or when the maximum depth of a tree is reached, one of the terminals is chosen.

2.4 Parameters of SR

Below, we list the complete start up parameters

Number of variables ( $n_{v a r}$ ) to use in terminals.
Number of generations ( $n_{g e n}$ ) in evolution of population.
Population size ( $n_{p o p}$ ) - At each generation, the population size is kept fixed.
Tournament size ( $n_{t o u r n}$ ) - A parent is chosen through a tournament from a randomly chosen subset of size $n_{t o u r n}$ from the population.
Flag whether to try double tournament $b_{d o u b l e t}$ - If $b_{d o u b l e t} = t r u e$ , each parent is chosen in separate tournaments. Otherwise, the best two in a single tournament are chosen as parents.
Maximum depth of tree ( $m_{d e p t h}$ ).
Copy reproduction ( $n_{c o p y}$ ) - best $n_{c o p y}$ individuals of a population are simply copied to the next generation. These individuals may still participate in sexual reproduction.
Fraction of new children produced through sexual reproduction in population ( $P_{x o v e r}$ ) - Number of children in each generation is $n_{p o p} \times P_{x o v e r}$ .
Mutation probability ( $P_{m u t}$ ) - Mutation probability per node
Drop probability ( $P_{d r o p}$ ) - Probability per node to drop it in a mutation
Probability to pick a constant ( $P_{c o n s t}$ )
Probability to pick $i$ th function or operator ( $P_{o p, i}$ )

For the examples in the next section, some of the parameters of interest are: $P_{x o v e r} = 0.5$ , $P_{m u t} = 0 0.1$ , $n_{p o p} = 5002500$ , $n_{t o u r n} = 0 0.2 \cdot n_{p o p}$ . The basic arithmetic operators $(+, -, \times, \div)$ have probabilities between 5% to 25% each. The other functions and operators take up the remaining probability equally. These probabilities vary for each run.

2.5 Initialization, training, fitness evaluation and evolution of SR

In the initialization stage, random expressions are created to populate the first generation. The user can choose to create an initial population with the expression filled to the maximum tree depth or not. This depends on the problem being considered. For the purpose of mass reconstruction, we find that having the initial population filled to the maximum depth of the tree is not necessarily beneficial, since most of the initial population has poor fitness and get discarded after a few generations.

No distinction is made between the training and testing sample. One half of the input sample is randomly chosen for the training sample and the other half is the testing sample. This choice is made for every generation. Fitness of each individual is evaluated on the testing sample according to the criteria for a problem.

Evolution of population is done through two main avenues. First is the sexual reproduction, where parents are not chosen randomly, but through a tournament. In a tournament, random subpool of $n_{t o u r n}$ individuals are chosen randomly and the best two individuals are chosen as parents in a single-tournament method. In a double tournament method, the parents are chosen from two tournaments by selecting the best fit individual in each tournament.

Once the parents are selected, then the subnodes are selected from two parents and swapped in place. We accept only one child from a set of parents, this choice is random. It is well-documented in literature that selecting a random subexpression for swapping leads to trivial or minimal modifications as nodes that have terminals are chosen more often. Therefore, we choose internal nodes 90% of the times for swapping.

Asexual reproduction is done by having mutations occurring in random nodes in the expression tree. Selected node is either replaced with another expression or it is dropped. If it is to be replaced, an expression tree randomly generated using the parameters of the symbolic regression algorithm. An individual which successfully had a child can still participate in an asexual reproduction. We finally, let the “elite” ( $n_{c o p y}$ ) to go to the next generation without modifications.

Figure 2: Creation of a child through cross-over operation. The cross over of genes occuring between $v 3$ and $v 4 - v 1$ from the two parents yields a new child $(v 1 + v 2) * (v 4 - v 1) .$

2.6 Evaluation of Fitness

Fitness function is problem specific. Evaluation of fitness is done on a set of test sample by traversing the expression tree. This is computationally the most time consuming part of SR. Expressions which do not yield machine-real numbers are assigned a number corresponding to a very poor fitness, marking them likely for removal. An individual with a more compact expression is favorable and a small penalty proportional to the number of nodes is imposed on each individual.

2.7 Dimensionally Constrained Symbolic Regression

Symbolic regression must yield a dimensionally sensible result if it is to be interpretable. In this application, we would like to construct equations related to mass of a particle. In natural units ( $c = 1$ and $ℏ = 1$ ), the mass has the same physical dimensions as momentum and energy. We will assume that the target dimension is $p^{n}$ , i.e., momentum raised to integer power, and that variables have dimension 1 or 0.

We see that the arguments to operators or functions must satisfy certain constraints. If $d (a)$ stands for the physical dimension of expression $a$ :

$a \pm b$ : $d (a \pm b) = d (a) = d (b)$ .
$a \times b$ : $d (a \times b) = d (a) + d (b)$ .
$a \div b$ : $d (a \div b) = d (a) - d (b)$ .
$\sqrt{a}$ : $\frac{d (a)}{2} = d (\sqrt{a})$
$x^{a}$ : $d (x^{a}) = d (x) * a$ and $d (a) = 0$
Transcendental function $f (a)$ : $d (f (a)) = 0$ and $a = 0$ .

Note that the argument to a transcendtal function must be dimensionless.

Using these rules, we can impose that an expression built by symbolic regression should have the correct physical dimensions. We would need to reformulate the transformation rules for evolution compared to the regular SR.

An $n$ (an integer) dimensional term can be built recursively by applying one of the following rules:

$n - 1$ -dim term $\times$ $1$ -dim term
$n$ -dim term $\pm$ $n$ -dim term
$n$ -dim term $\times$ $0$ -dim term
square root of an $2 n$ -dim term

A 0-dimensional term can be created by dividing two 1-dimensional terms or applying a transcendental function to argument of 0-dimension. And -1-dimensional term can be created by taking an invervse of 1-dimensional term. One can see that any integer dimensional term can be created using these rules.

In DCSR, rules for cross-over and mutation need to be modified accordingly. Exchange or replacement of expressions can occur only among those with the same physical dimensions, otherwise it is not permitted.

3 Application of Symbolic Regression to Mass Measurement Problem

In this section, we will apply symbolic regression to a few problems of mass determination in high-energy physics. In the first two cases, we apply it to the problem where the answer is well-known. And in the third case, we try to find a symbolic expression for the Higgs mass that can be applied to a case where Higgs boson decays to $W^{+} W^{-}$ and each $W$ boson decays leptonically.

3.1 Invariant mass

As a first example, we choose a trivial example, one for which we already know the answer. The variables given are the components of 4-vector,

(p_{x}, p_{y}, p_{z}, E)

, and the target value we would like to obtain is

m^{2} = E^{2} - p_{x}^{2} - p_{y}^{2} - p_{z}^{2}

. A 1000 toy data sample was generated satisfying the relationship.

Figure 3: Evolution of minimum absolute error of symbolic regression trained with invariant mass data. Out of 100 runs, 26% yield the correct formula $m^{2} = E^{2} - p_{x}^{2} - p_{y}^{2} - p_{z}^{2}$ . Evolution of fitness of 36th run is shown by itself.

Symbolic regressions with different values of parameters were tried on this data. Fitness function is $\frac{1}{n} \sum_{n} | m - m_{t r u e} |$ , where $n$ is the size of data. Figure 3 shows evolution of fitnesses of the best fit individuals in each run. In 26% of the runs, the correct formula for invariant mass is found. The right plot of Fig. 3 shows the evolution of the best fit individual of run 36. At the start of the run, the best fit expression is $\sqrt{| E^{2} cos p_{z} - E |}$ , which is meaningless physically. By generation 86, the equation $m^{2} = E^{2} - p_{x}^{2} - p_{y}^{2} - p_{z}^{2}$ is found. When DCSR is used on the same sample, the chances for success increases to 68%.

3.2 Events with Missing Information - Transverse Mass

A more realistic and interesting application of symbolic regression is to a problem where information is lost. A classic example would be where a massive particle decays into a charged lepton and a neutrino in a hadronic collision. Transverse mass $M_{T} = \sqrt{p_{T ℓ} ⧸ ⧸ E_{T} - {\to p}_{T ℓ} \cdot ⧸ {\to p}_{T}}$ is typically used to measure the mass of the $W$ boson. To test whether we can extract this mass relationship, a phase-space generator is used to create the momenta of lepton and neutrino from a hypothetical massive particle with an exponentially falling initial transverse momentum distribution for particle masses between 0 GeV and 100 GeV.

Figure 4: Evolution of minimum of fitness in 100 runs without DCSR (left) and with DCSR (right) during training using lepton and neutrino momenta variable set as inputs $(p_{T}, p_{x}, p_{y}, p_{z}, ⧸ E, ⧸ E_{T x}, ⧸ E_{T y})$ . With DCSR, the fraction of runs that find the equation for transverse mass increases dramatically.

The fitness function used was $\frac{1}{n} \sum^{n} | m_{i, e s t} - m_{i, t r u e} | / m_{i, t r u e}$ , where $m_{i, e s t}$ is the value returned by $i^{t h}$ individual. Some trees reach a global minimum as shown in Fig 4. For those that reach a global minimum, the equation is $m_{e s t}^{2} = 2.38 (p_{T ℓ} ⧸ ⧸ E_{T} - {\to p}_{T ℓ} \cdot ⧸ {\to p}_{T})$ , which is the same as the definition of $M_{T}$ except for the overall constant factor.

Depending on the parameters of the genetic algorithm, it may take longer to reach this minimum. Some may fail to reach this minimum For example, if the minimization criteria were

\sqrt{\frac{1}{n} \sum^{n} (m_{i, e s t} - m_{i, t r u e})^{2}}

, almost all trials of symbolic regression fail to reach the global minimum. The populations get quickly trapped into a local minima. This is the result of using a minimizing criteria that heavily penalizes the outliers. It is a known feature of genetic algorithm that having individuals with poor fitness in the gene pool is important in population to reach the global minimum. Without DCSR, 7% of the trials find the solutions, while with DCSR the success rate increases to 40%. We find that constraining the terms allowed increases the speed of convergence and the chances for a successful convergence.

3.3 Mass Sensitive Variable in $H \to W W^{*} \to ℓ^{+} ν ℓ^{-} ¯ ν$

Higgs mass determination in $H \to W W * \to ℓ^{+} ν ℓ^{-} ¯ ν$ in hadron colliders is an inportant problem. In this channel, two lepton momenta ${\to p}_{ℓ 1}, {\to p}_{ℓ 2}$ and the vector sum of the two neutrino transverse momenta $⧸ {\to E}_{T} = ⧸ {\to E}_{T ν 1} + ⧸ {\to E}_{T ν 2}$ are measured in experiments. Since there are only two equations related to neutrino momenta, the system is under-constrained. If we knew both $W$ bosons were real, we would still need two extra equations to constrain the system. Therefore one cannot solve for the neutrino momenta exactly even in principle.

Existing studies relied on analysis of kinematics to find expressions that behave linearly to the Higgs boson mass [3, 5]. In this study, we apply symbolic regression to the problem and find an expression that not only shows linear behavior, but also whose widths of the mass distribution are narrow.

Symbolic regression is applied to a data generated with PYTHIA $p p \to H \to W W^{*} \to ℓ^{+} ν ℓ^{-} ¯ ν$ at $\sqrt{s} = 14$ TeV with $M_{H}$ varying from 120 GeV to 200 GeV [8]. Detector simulation is not performed on this data in order to find the ideal expression. Momentum components and energy of the two charged leptons ( $p_{1 T}, p_{1 x}, p_{1 y}, p_{1 z}, E_{1}, p_{2 T}, p_{2 x}, p_{2 y}, p_{2 z}, E_{2}$ ) and missing $E_{T}$ information ( $⧸ E_{T}, ⧸ E_{x}, ⧸ E_{y}$ ) are used as input variables for the symbolic regression. The fitness function used is the average of fractional absolute difference: $\frac{1}{N} \sum_{i} | M_{r e c, i} - M_{H, i} | / M_{H, i}$ .

Without DCSR, the symbolic regression is not able to yield meaningful results. This seems to be due to the greater number of variables used. The number of terms of dimension 2 with only multiplication allowed is 78, which makes the possible function space to explore very large.

If fractional root mean-squared (RMS) ( $\frac{1}{N} [\sum_{i} (M_{r e c, i} - M_{H, i})^{2} / M_{H, i}^{2}]^{1 / 2}$ ) were used as the fitness function, the symbolic regression would get trapped into local minima even with DCSR since outliers pay a heavy penalty.

Figure 5

shows evolution of fitness of best-fit individuals in 100 runs as a function of the number of generations. DCSR is able to converge on meaningful results and yields the best estimate for the

M_{H}^{2}

S_{m a s s}^{2} = 2 p_{1 T}^{2} + 2 p_{2 T}^{2} + 3 (p_{1 T} p_{2 T} + ⧸ E_{T} (p_{1 T} + p_{2 T}) - ⧸ {\to E}_{T} \cdot ({\to p}_{1 T} + {\to p}_{2 T}) - 2 {\to p}_{1 T} \cdot {\to p}_{2 T}) .

Symmetry of the two leptons in the system is recognized by the symbolic regression automatically, even though symmetry condition was not imposed. More detail on the Higgs mass measurement using this result is discussed in Ref. [9].

Figure 5: Left: Evolution of best fit individuals function in 100 runs. Right: Distribution of ratio of predicted mass to true mass, $m_{p r e d} / m_{H}$ , versus the true Higgs mass $m_{H}$ . The sample was generated using PYTHIA.

3.4 Computational performance

With our implementation, on a 24 thread X5650 Xeon machine, it takes approximately 1 day to complete the 100 runs for the Higgs mass example. Most of the time is spent in evaluation of the fitness function. We noticed memory fragmentation occurring and evaluation tended to become slower as time went on. A Mathematica implementation, which was simple to write, did not suffer from memory problems. However, it was inherently slow. Using the parallel processing capabability of Mathematica 7, even when we used all 24 threads, one run took about a day to complete.

With the advent of higher multiplicity core CPUs in the future, it may become computationally favorable to use this tool. One avenue that needs to be explored is to consider the backgrounds when optimizing. At the moment, it would be very expensive to evaluate fitness including backgrounds since a simple fitness function would no longer suffice.

4 Conclusion

We have defined and implemented dimensionally constrained symbolic regression for use in high-energy physics problems. We find that it reproduces some of the well-known results, such as invariant mass and transverse mass. We used it to construct a function that is linear and sensitive to the Higgs boson mass $M_{H}$ . that has not been previously known. We expect that this method would be useful in cases with more than one undetected particles, where it is not trivial what variable could give the best performance.

References

[1] G. Altarelli et al.Nucl. Phys. B 246 (1984) 12.
[2] J. R. Koza, 1992. Genetic Programming: On the Programming of Computers by Means of Natural Selection, Cambridge, MA: MIT Press.
[3] Alan J. Barr, Ben Gripaios, Christopher Gorham Lester, J. High Energy Phys. 0907 (2009) 072.
[4] M. Dittmar, H. K. Dreiner, Phys. Rev. D 55 (1997) 167.
[5] Kiwoon Choi, Suyong Choi, Jae Sik Lee, Chan Beom Park, Phys. Rev. D 80 (2009) 073010.
[6] Holland, John H. 1975. Adaptation in Natural and Artificial Systems, University of Michigan Press.
[7] http://root.cern.ch
[8] Torbjörn Sjöstrand, Stephen Mrenna and Peter Skands, J. High Energy Phys. 0605 (2006) 026.
[9] S. Choi, 1006.4998.
[10] Wolfram Research, Inc., Mathematica, Version 7.0, Champaign, IL (2010).

Dimensionally Constrained Symbolic Regression

Authors

AI Feynman: a Physics-Inspired Method for Symbolic Regression

Globally Optimal Symbolic Regression

Discrete symbolic optimization and Boltzmann sampling by continuous neural dynamics: Gradient Symbolic Computation

A general methodology to assess symbolic regression algorithms using the generation of random equations with uniform random sampling

Where are we now? A large benchmark study of recent symbolic regression methods

Symbolic Regression Driven by Training Data and Prior Knowledge

Does Symbolic Knowledge Prevent Adversarial Fooling?

1 Introduction

2 Dimensionally Constrained Symbolic Regression and Its Implementation