Deep Learning Analysis of Deeply Virtual Exclusive Photoproduction

The cross section for $ep\rightarrow e^{\prime}p^{\prime}\gamma$ reads,

where $T$ is a superposition of the DVCS amplitude and it’s background process (Bethe-Heitler) BH,

The cross section is a differential form of five independent variables, the four momentum transfer squared between the initial ($e$) and final ($e^{\prime}$) electrons, $Q^{2}$; $x_{Bj}=Q^{2}/(2M\nu)$ where $\nu$ is the energy transfer between the initial and final electrons and $M$ is the proton mass; $t$, the (negative) four momentum transfer squared between the initial and final proton; the angle $\varphi$ between the lepton and hadron planes, and $\varphi_{S}$, the azimuthal angle of the transverse proton spin vector. We will consider in our analysis available fixed target data for an unpolarized target where the measured set of observables includes

A Deep Neural Network (DNN) is an effective function approximation method consisting of several layers of computation, in which each layer computes a new representation of the previous layer’s outputs. In the simplest case, each layer consists of an affine transformation of the previous layer’s representation of the data, followed by an element-wise nonlinearity. For example, the output of the first layer is:

Where $z$ is a nonlinear activation function applied to each element of its input vector. $\mathbf{W}$ is an $m\times n$ matrix where $m$ is the dimension of the previous layer’s representation and $n$ is the dimension of the current layer’s output. $\mathbf{b}$ is a vector of size $n$. The entire model can be written as a composition of functions

Where $y$ is the true value of our predicted variable, often called the ‘label’ or ‘target’ value. In our case, backpropagation is used to compute the partial derivatives of the loss function on a batch of $B$ samples w.r.t each of the parameters in $\theta$, and the resulting gradient vector is used to update the parameters via Stochastic Gradient Descent (SGD):

Where $\lambda$ is a learning rate parameter. Modern optimizers (Kingma and Ba, 2014) (Hinton et al., ) often dynamically adjust $\lambda$ based on heuristics like momentum, where consecutive gradient steps in a similar direction increase the learning rate and rapid changes in direction push it towards zero.

We compiled a dataset from several DVCS experiments Georges (2018); Defurne et al. (2015, 2017); Jo (2012), covering the range of kinematic bins shown in Figure 1. We discarded data sets that listed no explicit $\varphi$ dependence. We used a common format for collecting and recording our data sets such that they list the kinematic dependence of the cross section points, the phi dependence, the polarization of the observable, and the statistical and systematic errors.

Each entry consists of the input kinematic variables ($x_{Bj}$, $t$, $Q^{2}$, $E_{b}$, $\varphi$) as well as the value and experimental error of the cross section. In total, we are working with 3,862 unpolarized and 3,884 polarized data points. Because these measurements are typically taken by sweeping across $\varphi$ values in a fixed kinematic bin, most points are unique only in $\varphi$, meaning there are relatively few values of $x_{Bj},t,Q^{2},E_{b}$. This presents an interesting challenge for ML approaches, which become much more capable of predicting across $\varphi$ values than they are of extrapolating these patterns to distant kinematic space.

To help measure this effect, we use a $70\%/20\%/10\%$ split to separate the kinematic bins into three groups, which we’ll call the training set, validation set and test set, respectively. Note that this differs from the standard ML practice of splitting the datapoints themselves into random groups. The training set is used to directly adjust the parameters of the model. The validation set is used as an indicator of performance during training - either for early stopping (see Sec III.3), or for hyperparameter selection¹¹1A hyperparameter is a fixed constant that controls some aspect of an algorithm’s behavior. They can have a significant impact on performance, but cannot be learned like many of the model’s core parameters. As an example, SGD optimizes the parameters of a model, but the initial learning rate and batch size of the optimizer are predetermined hyperparameters.. Final performance will be reported on the test set, which is made up of kinematic bins with which the model has never been trained or evaluated. A model that performs well on the test set is capable of predicting entire $\varphi$ curves in unseen regions of kinematic space, under the assumption that we are testing on inputs drawn from the same distribution as our training data. This assumption often breaks down in practice, as we’ll see in Section III.4. Generalizing to unseen kinematic bins is a more difficult challenge than generalizing to held-out datapoints in bins that were part of the training set, because the overall shape of the $\varphi$ curve can be quite predictable as long as nearby points can provide a good frame of reference.

Based on the analysis of (Kriesten et al., 2020), we experimented with transforming the dataset of ($x_{Bj}$, $t$, $Q^{2}$, $E_{b}$, $\varphi$) by concatenating (sin($\varphi$), cos($\varphi$), sin($2\varphi$), cos($2\varphi$), sin($3\varphi$), cos($3\varphi)$). These harmonics could serve as additional features that assist the model in learning the cross section’s $\varphi$ dependence. We’ll refer to the use of this additional input information as “Harmonic Features.”

We also scale the cross section values by 1000 to mitigate floating-point underflow and stabilize training. Unless otherwise noted, the plots below have been re-scaled to their original units.

First, we examine the possibility of applying basic regression techniques to cross section fits to get an understanding of the difficulty of the problem and establish a performance baseline. We use the implementations provided in Pedregosa et al. (2011):

• $\bullet$

  Linear is a simple linear regression model. A vector of weights $\mathbf{W}$ and a scalar bias $b$ map the input kinematics, $x$, to the cross section prediction with $\hat{\sigma}=\mathbf{W}x+b$. These parameters are the solution to a closed-form system of equations that minimizes the mean squared error (MSE) of the model’s predictions Kenney and Keeping (1962).

• $\bullet$

  SVR is a Support Vector Machine Regression model Vapnik (1999). The model parameters, $\mathbf{W}$ and $b$ are the solution to the optimization problem Chang and Lin (2011):

  	$\displaystyle\mathop{min}_{\mathbf{w},b,\gamma,\gamma^{*}}$	$\displaystyle\frac{1}{2}\mathbf{w}^{T}\mathbf{w}+C\sum_{i=1}^{l}\gamma_{i}+C\sum_{i=1}^{l}\gamma^{*}_{i}$
  	subject to	$\displaystyle\mathbf{w}^{T}\rho(\mathbf{x}_{i})+b-\sigma_{i}\leq\epsilon+\gamma_{i},$
  		$\displaystyle\sigma_{i}-\mathbf{w}^{T}\rho(\mathbf{x}_{i})-b\leq\epsilon+\gamma^{*},$
  		$\displaystyle\gamma_{i},\gamma^{*}_{i}\geq 0,i=1,\dots,l.$

  Where $C$ and $\epsilon$ are hyperparameters that control error tolerance and $\rho$ is a kernel function that is used to map input vectors to a space where they may be more linearly separable. We experiment with both polynomial and radial basis function kernels.

We train two separate NNs, one to predict $\sigma_{UU}$ and the other to predict $\sigma_{LU}$. There are several major challenges in applying Deep Learning to this problem. The first is the range of the cross section data values; the target values for our networks’ predictions span seven orders of magnitude in both the polarized and unpolarized cross sections. While the last layer of our architecture uses a linear activation function that can output any scalar, this wide range in targets destabilizes training. When kinematic points with large cross sections are sampled in a batch of training data, they can create outlier loss values that increase the size of our update gradients and make large adjustments to model parameters that decrease performance in regions with a smaller cross section. This effect is magnified by the MSE function, which is quadratic in the absolute error of our models’ predictions. We address this by changing our loss function to the Huber loss:

The Huber loss is quadratic in the model’s error when that error is less than a hyperparemter $\delta$, and linear otherwise. This limits the effect of large cross section values relative to MSE. The median value of the cross section data is relatively small, so the squared error term is stable for most datapoints and therefore the specific choice of $\delta$ is unimportant; we fix $\delta=1$. Many data points have error bars that are listed well above $100\%$ of the cross section value. Fitting to the center of such a wide range can create irregularities that generalize poorly and encourage overfitting. Ideally, our model would pay less attention to highly uncertain measurements. Towards this goal, we consider a variation on the standard Huber Loss that incorporates the relative size of the experimental error bars for each datapoint. Let $\Delta\sigma$ represent the size of the symmetric error bar, while $\Delta\sigma_{\downarrow}$ and $\Delta\sigma_{\uparrow}$ represent the negative and positive asymmetric error, respectively. The Error Weighted Huber Loss (EW Huber) is then defined:

The EW Huber discounts the loss value of datapoints where the total size of the error bar is large relative to the scale of the cross section. Points with no experimental error stay unchanged and the effect of the weighting increases with the experimental uncertainty.

In many regions of kinematic space, the unpolarized cross section is dominated by the Bethe-Heitler component, which can be computed as a known function of the network’s inputs. We take advantage of this by implementing the Bethe-Heitler calculation as a custom neural network layer with fixed parameters that cannot be adjusted by gradient descent. This layer computes the Bethe-Heitler component of the input point, which is then multiplied by the output of the standard NN architecture to produce the predicted cross section. This effectively changes the function the NN must learn to approximate from $\sigma$ to the simpler $\sigma/BH$. We test these approaches on increasingly large NN sizes to see how they affect performance, and the results are shown in Figure 4. The Huber loss turns out to be essential for successful training. Performance improves with model size - perhaps the biggest advantage of Deep Learning techniques. The BH Ratio simplifies the NN function enough to help small models, but as models grow larger their representational capacity improves and this simplification becomes unnecessary.

While the cross section measurments used in our analysis consist of around 4,000 datapoints, that total will likely see significant growth from future experiments. We are interested in developing a method that can scale as new measurements are collected. This can be tested by expanding the kinematic region to include a larger range of $x_{Bj},t$, and $Q^{2}$. In particular, we generate pseudo-data in regions outside the range where data currently exists to measure the network’s compatibility with theoretical calculations of the cross section.

When making predictions that have an effect on decisions or later calculations, we’d like to know the confidence level of our model - how sure are we that our predictions are accurate? This confidence is quantified as model uncertainty. Recent literature (see Refs.Kumericki et al. (2011); Moutarde et al. (2019) for DVCS, Ball et al. (2017, 2012) and references therein for inclusive scattering) has adopted the “replica method” as a means of estimating uncertainty. The replica method is similar to the Bootstrap aggregating (Breiman, 1996) approach that is common in ML, in which a set of models is trained by sampling from the dataset with replacement and using the disagreement of their predictions to estimate uncertainty. The replica method consists of running the learning algorithm many times on versions of the dataset where points are randomly perturbed relative to their experimental error bars. This is meant to propagate the experimental error into the results of the model. After many cycles of resampling and retraining, we have a set of models trained on slightly different datasets. The mean output of this set of models is used as the prediction, and the standard deviation serves as an approximation of uncertainty. While this is a valid way to incorporate experimental/epistemic uncertainty into the modeling process, the repeated training can be prohibitively expensive when working with large models or datasets. Using multiple copies of the model parameters can be memory-intensive, difficult to vectorize, and less computationally efficient when making predictions. Ideally, we could recover model uncertainty from a single neural network.

One possibility is a Bayesian Neural Network, where we optimize a distribution over weights and sample from that distribution between each prediction Blundell et al. (2015). While effective, these networks have significantly more parameters and can be more difficult to train. A simple alternative is to leave dropout turned on during inference. By resampling our binary mask over the weights between each forward pass, we create a prediction distribution that can be used to approximate model uncertainty. This technique can be shown Gal and Ghahramani (2015) to be an approximation of a Gaussian Process. Our model generates a batch of $k$ predictions over the same datapoint, each with different dropout masks. Unless otherwise noted, we set $k=1000$. This smooths the uncertainty bands for plotting, but this value can be decreased if inference speed is a priority.

First, we examine the model’s ability to fit the range of $\varphi$ values in existing kinematic bins ($x_{Bj}$, $t$, $Q^{2}$, $E_{b}$). FemtoNet’s predictions in unseen kinematic bins from the test set are shown in Figure 10. The models learn to predict the cross section within experimental error bars, and are robust to noise and outliers, particularly in the Hall B datasets.

The total unpolarized cross section at leading twist contains 3 pieces as shown in Eq.(2): a piece that is bilinear in Compton form factors (DVCS), an interference piece that is linear in Compton form factors ($\mathcal{I}$), and a piece that has no dependence on Compton form factors and is an exactly calculated background parametrized by the elastic form factors (BH). The Bethe-Heitler cross section term, in most kinematic regimes, dominates the total cross section; often an order of magnitude or larger than the interference and DVCS terms combined. We can calculate the background to a well known degree of accuracy using covariant four vector products. The error and the noise of the cross section is thus due entirely to the terms which are parametrized by the Compton form factors. Since the background dominates the cross section, it then also dominates the noise inherent in scattering cross section data, effectively reducing the errors and simplifying the functional form. It is demonstrated in Table 2 that the neural network can predict the unpolarized cross section with a relatively small median percentage discrepancy.

This is in contrast to the polarized cross section which has a drastically different result. The functional form of the polarized cross section we find is more difficult for the neural network to learn. We can explain this by the presence of a large contribution of noise in the polarized data. By parity constraints, the longitudinally polarized beam on an unpolarized target (LU) cross section has no Bethe-Heitler background nor a term bilinear in Compton form factors at leading twist. It is linear in Compton form factors; as such, there is no error / noise suppression as in the unpolarized case.

We can compare the LU polarized cross section to an equivalent reduced unpolarized cross section where we have subtracted off the Bethe-Heitler background to investigate the impact of the dominant term on the errors as shown in Figure 11. Notice that the unpolarized cross section tends to converge to a relative error similar to that of the polarized cross section at similar kinematics. Where the relative error is defined as

The background in the unpolarized cross section tends to simplify the functional form of the cross section while simultaneously smoothing out the noise and decreasing relative error. This, however, presents a problem for extracting Compton form factors where the cross section terms they parametrize tend to have large noise components.

We aren’t restricted to kinematic bins that have already been measured - the model learns to extrapolate between existing bins. We can sweep through regions where data does not exist, and use the model’s uncertainty estimate to advise our use of its predictions. In Figure 12 we illustrate this point for two particular kinematic values.

This result is consistent with what reported in the upper panel of Fig.7 where the model’s uncertainty in the experimental data domain is very low. Notice that the neural network seems to be able to predict the $t$ dependence much more accurately than the $xi$ dependence. This result could mean, on one side, that the $t$ dependence has a simpler and consistent functional form throughout the phase space than the $\xi$ dependence. The latter will, in fact, be affected by perturbative QCD evolution. On the other hand, the kinematic spread of the data for the $t$ variable covers a wider section than for $\xi$: this feature might also influence the outcome. Clearly, the availability of more data in a wider range in $\xi$ for fixed $t$ will help interpreting this result.

While the dataset contains thousands of datapoints, it contains significantly fewer unique kinematic bins; most of the individual points come from kinematic sweeps in the azimuthal angle $\varphi$. The impact of this can be seen empirically in the model’s uncertainty estimates. The uncertainty grows quickly as we move between different kinematic regions, while it remained relatively constant between different $\varphi$ values in a fixed region.