Artificial neural network modelling of generalised parton distributions

We discuss the use of machine learning techniques in effectively nonparametric modelling of generalised parton distributions (GPDs) in view of their future extraction from experimental data. Current parameterisations of GPDs suffer from model dependency that lessens their impact on phenomenology and brings unknown systematics to the estimation of quantities like Mellin moments. The new strategy presented in this study allows to describe GPDs in a way fulfilling theory-driven constraints, keeping model dependency to a minimum. Getting a better grip on the control of systematic effects, our work will help the GPD phenomenology to achieve its maturity in the precision era commenced by the new generation of experiments.


Introduction
Generalised parton distributions (GPDs) [1][2][3][4][5] are widely recognised as one of the key objects to explore the structure of hadrons.They encompass information coming from one-dimensional parton distribution functions (PDFs) and elastic form factors (EFFs).GPDs allow for a hadron tomography [6][7][8], where densities of partons carrying a fraction of hadron momentum are studied in the plane perpendicular to the hadron's direction of motion.GPDs also provide access to the matrix elements of the energy-momentum tensor [2,3,9], making it possible to evaluate the total angular momentum and "mechanical" properties of hadrons, like pressure and shear stress at a given point of space [10][11][12][13].
Experimental access to GPDs is mainly possible thanks to exclusive processes occurring on hadrons remaining coherent after a hard scale interaction.Some notable processes of this type are deeply virtual Compton scattering (DVCS) [3], time-like Compton scattering (TCS) [14] and deeply virtual meson production (DVMP) [4,15].All of them allow to study transitions of hadrons from one state to another, with a unique insight into changes taking place at the partonic level.GPDs have been primarily studied in leptoproduction experiments, in particular those conducted in JLab, DESY and CERN, and are key objects of interest in programmes of future electron-ion colliders, like EIC [16,17], EicC [18] and LHeC [19].Many data sets have already been collected and reviewed for DVCS [20,21] and DVMP [22], while the first measurement for TCS has been completed recently [23].
Although several datasets are already available for fits, the extraction of GPDs is far from being satisfactory.The main reasons are: i ) sparsity of available information.GPDs are multidimensional functions, so one needs much more data to constrain them in the phase-space of kinematic variables, comparing to e.g.one-dimensional PDFs.Furthermore, one needs to cover this phase-space by data collected for various processes and experimental setups, which is required to distinguish between many types of GPDs and contributions coming from various quark flavours and gluons.ii ) complexity of extraction.In order to fully benefit from available sources of information about GPDs, such as data collected in exclusive measurements, PDFs and EFFs for boundary conditions, lattice QCD, one needs to know and implement links connecting those sources with GPDs.The extraction of GPDs requires a deconvolution of amplitudes measured in exclusive processes, which is not trivial and in some cases does not even possess a unique solution [24].Both the evolution of GPDs and the description of exclusive processes must be understood to a level allowing for a robust extraction.This requires a substantial effort and a careful design of tools aggregating GPD-related theory developments.iii ) model dependency.Currently available phenomenological GPD models, like GK [25][26][27] and VGG [28][29][30][31], use similar Ansätze with a rigid form and therefore cannot be considered as diverse sources for the estimation of model uncertainty.The severity of such a model dependency and its impact on e.g. the extraction of orbital angular momentum have never been studied in a systematic way, partially because there were no tools to do so.
The answer of the GPD community to those problems so far has been as follows.More data sensitive to GPDs will be delivered by a next generation of experiments.These experiments, performed in many laboratories and with multiple setups, will provide the much needed input to the GPD phenomenology.We witness a substantial progress in the description of exclusive processes and in the systematic use of GPD evolution equations [32].The development of lattice QCD techniques is accelerating, as proved by Ref. [33].The PARTONS [34] and GeParD [35] projects provide aggregation points for GPD-related developments and practical know-how in phenomenology methods.Techniques to stress model dependency are developed for the extraction of amplitudes, the latter used to access the "mechanical" properties of nucleons [36,37].However, to the best of our knowledge, no substantial progress has been recently made to address the problem of model dependency of GPDs extractions.
In this Article, we discuss new ways of modelling GPDs which provide flexible parameterisations inspired by artificial neural network (ANN) techniques.These parameterisations may be used as tools to study the model dependency affecting extractions of GPDs and derived quantities, like orbital angular momentum of hadrons.We address aspects of a practical nature, e.g.dealing with many free parameters.Our work aims in providing a better grip on the control of systematic effects, much needed in front of the precision era of GPD extractions.
The Article is organised as follows.In Sect. 2 we provide the theoretical background for our line of research.The modelling of GPDs in (x, ξ)-space is given in Sect.3, while that in (β, α)-space is discussed in Sect. 4. We provide the summary in Sect. 5.

Theoretical background
GPDs are functions of three variables: the average longitudinal momentum fraction carried by the active parton, x, the skewness ξ, which describes the longitudinal momentum transfer, and the four-momentum transfer to the hadron target, t.GPDs also depend on the factorisation scale, µ F , being the variable entering evolution equations.A knowledge of these equations allows to run the evolution starting from a reference scale, µ F,0 , where GPD models are defined, to any µ F .We therefore only consider GPDs at µ F,0 , and in the following we suppress the µ F dependence for brevity.For illustration we will only consider the GPD H(x, ξ, t) for a quark of unspecified flavour, but we note that the discussion can be easily extended for other GPDs.H(x, ξ, t) is a real and ξ-even function as a consequence of QCD invariance under discrete symmetries.Without loss of generality, we only consider positive ξ in this work.
GPD models must fulfil a set of theory driven constraints: i ) reduction to the PDF q(x) in the forward limit: ii ) polynomiality [38][39][40], which is required to keep GPDs invariant under Lorentz transformation.The property states that any Mellin moment of a GPD, A n , is a fixed-order polynomial in ξ: where we call A n,j (t) a Mellin coefficient.We note that A 0 (ξ, t) = A 0,0 (t) ≡ F 1 (t) is the Dirac EFF form factor.The polynomiality "entangles" the x and ξ dependencies.iii ) positivity constraints [41][42][43][44][45][46][47][48][49], which are inequalities ensuring positive norms in the Hilbert space of states.Because of the variety of these inequalities, we do not quote them all here.Instead, we refer the Reader to the aforementioned references, and only note that the constraints typically involve several types of GPDs.In the following we will illustrate how to deal with the positivity constraints with this simple inequality [41,42], which is applicable in the DGLAP region, i.e. for x > ξ.
For the sake of completeness, we write a direct relation between Mellin moments, which are extensively discussed in this Article, and conformal moments given by: where C 3 /2 n (x) are Gegenbauer polynomials: with coefficients A finite number of Mellin moments is needed to express a single conformal moment of a fixed order, and vice versa, a finite number of conformal moments is needed to express a single Mellin moment, The conformal moments conveniently diagonalise the evolution equations at leading-order (LO) [50], but they also appear in the modelling of GPDs, like for instance in Ref. [51].
Working in (β, α)-space allows us to relatively easily fulfil both the reduction to PDFs and the polynomiality property.However, any projection of experimental observables using GPD models defined in that space, or vice versa, any attempt to find free parameters of such models from experimental data, requires either the Radon transform or its inverse counterpart.This is typically done numerically, which severely slows down computations.In some cases the Radon transform can be performed analytically, but it requires a rather simple form of double distributions and may be numerically unstable due to delicate cancellations of small numbers, typically occurring at small ξ (see e.g. the O ξ −5 factor in Eq. ( 27) of Ref. [53]).The inverse Radon transform poses a challenge of its own, as reported for instance in Ref. [54].
We close this section with the remark that the core of this work is related to preserving proper forward limits, polynomiality and positivity of GPDs models.This primarily does not involve the t variable, which in general is not relevant in our study and as in other analyses of this type will suppressed for brevity.In summary, we are left with the problem of modeling GPDs in 2D spaces: either (x, ξ) or (β, α).

Principles of modelling
The polynomiality property suggests that the A n,j coefficients introduced in Eq. ( 2) provide a convenient basis to describe GPDs in a flexible modelling based on multiple parameters.The coefficients are not expected to be correlated, making them true degrees of freedom for GPDs.Another reason to choose this basis is the relation between Mellin and conformal moments indicated in Sect. 2. This relation allows one to perform the LO evolution of GPDs in a fast and straightforward way.Finally, the strategy advocated here is a convenient way to incorporate calculations of Mellin moments performed by lattice QCD in the modelling of GPDs.
Because there are infinitely many A n,j coefficients, for a practical use we only consider a subset of them.This is justified by the fact that the evaluation of higher order Mellin moments becomes increasingly sensitive to the numerical noise: if H(x, ξ) is a bounded summable function, its Mellin moment A n (ξ) tends to zero when n gets to infinity, in which case it starts to drown in the noise in any evaluation.This prevents the extraction of high order Mellin moments from experimental data without adding prior assumptions.Evaluation on the lattice also suffers from increased uncertainty at higher orders.Therefore, even if in principle the full series of Mellin moments is required to unambiguously reconstruct a function, this situation almost never appears in practical situations.The general case to consider is actually when only a finite number of Mellin moments are known.An extra regularization is thus necessary to restore the uniqueness of the function defined by its Mellin moments (in the present case, this extra regularization is introduced by the choice of a polynomial basis or an ANN basis as described below).By linearity, we point out that this argument still holds true when dealing with conformal moments instead of Mellin moments.
In our modelling, we only let free the coefficients A n,j appearing in Mellin moments up to and including the order n ≤ N .The value of N is arbitrary and it controls the flexibility of models, just like the order of polynomials or the size of artificial neural networks used in modern parameterisations of PDFs.In other words, N fixes the number of degrees of freedom.For even values of N , which for brevity we only consider in this study, this number is (1 The coefficients A n,j for n > N and j ≤ N are recovered from those for n ≤ N by a specific reconstruction procedure of the x-dependence that is discussed in the following.On the contrary, we assume that the coefficients A n,j vanish for n > N and j > N .Note that this is true at µ F,0 where we are defining the model.For µ F = µ F,0 , evolution will make these once null coefficients to have non-zero values, but only for j ≤ n + mod(n, 2) since evolution preserves the polynomiality property.We summarise this paragraph with a sketch shown in Fig. 1.Fig. 1 The sketch indicates three groups of the polynomiality expansion coefficients, A n,j , appearing in the modelling after fixing the truncation parameter to N = 4.The coefficients that we assume to be constrained by data and lattice QCD are denoted by , those that are null, but only at the initial factorisation scale are denoted by , while those calculable from after fixing the recovery of x-dependence are denoted by .
The reconstruction of the x-dependence from a number of A n,j coefficients is the moment problem known in mathematics.Let us address it by expressing GPDs in a way suggested by the polynomiality property: Here, f j (x) is a function of x that satisfies the following relation to Mellin coefficients: where A n,j = 0 for j > n + mod(n, 2).Since H(x = ±1, ξ) = 0 and considering the polynomial form of Eq. ( 11), f j (x) has to vanish at the end points: Equations ( 12) and ( 13) are essential to constrain f j (x) and therefore to obtain valid GPD models.In the following we will try two bases for f j (x), i.e. two forms of this function: one based on monomials, which has already been studied in the literature (see Ref. [52] and references therein), and one inspired by the ANN technique.A given basis fixes the bridge between the picture of GPDs given by Mellin moments (ξ-dependent integrals over x) and the picture where GPDs are explicitly given as functions of x and ξ.We note that the selection of a given basis and the necessity of keeping a finite N introduces some model-dependency.We however expect this bias to become arbitrarily small for sufficiently large N .
It is important to note that GPDs may exhibit a singularity at x = ξ = 0, which corresponds to the singularity of PDF at x = 0.However, this singularity does not affect the evaluation of Mellin moments.Since here we are modeling GPDs with a polynomial in ξ, see Eq. ( 11), we are not able to effectively build a model that exhibits the singularity at only x = ξ = 0 (and not at x = 0, ξ = 0).It is a general flaw of this type of modelling, which motivated us to find ANNinspired GPDs in (β, α)-space, which we present in Sect. 4. Still, a direct modelling of (x, ξ)-space can be used for a subset of non-singular GPD models, like the pion model [55,56] used for a demonstration in Sect.3.4.In addition, one can consider a reference GPD model, H 0 (x, ξ), with a correct H 0 (x, ξ)/q(x) behaviour for ξ → 0, and constrain a polynomial model, H 1 (x, ξ), such that H 0 (x, ξ)+H 1 (x, ξ) correctly reproduces ξ > 0 data and e.g.lattice results for Mellin moments.This kind of hybrid modelling may be plausible for phenomenology, as models defined in (x, ξ)-space can be directly compared to data, and therefore their usage can severely speed up minimization procedures.

Polynomial basis
According to the Stone-Weierstrass theorem, any continuous function defined on a compact set can be uniformly approximated by a polynomial to any desired degree of precision.Applying this theorem to the twodimensional function H(x, ξ) suggests to express f j (x) that is defined in the |x| ≤ 1 interval by Here, w i,j are coefficients multiplying the x i ξ j monomials in the global expression of H(x, ξ).The order of polynomial ( 14) is N +2, and it is the minimal order required to satisfy both Eqs. ( 12) and ( 13).We note that the order may be higher, which provides extra degrees of freedom available in the modelling.Mellin coefficients can be evaluated from monomial ones in the following way: where we remind that N is even, and where Eq. ( 13) yields: To find the opposite relations we need to solve Eqs. ( 15), ( 16) and ( 17) for w i,j .We express all these equations in the matrix form: Here, A j = (A 0,j , . . ., A n,j , . . ., A N,j , 0, 0) T and w j = (w 0,j , . . ., w i,j , . . ., w N +2,j ) T , and the length of each vector is N + 3. The matrix C is: where: This matrix is only characterized by the parameter N .Therefore, as soon as this parameter is fixed the matrix is ready for the inversion (for instance numerically), and this can be done only once in the modelling procedure.We note that the presented realisation of GPDs in the monomial basis of minimal order N + 2 naturally arises in the context of the so-called dual parameterisation [10] after its initial formal representation has been recast as a convergent series.This parameterisation is physically motivated by an infinite series of tchannel exchanges.We also note that any attempt of modelling GPDs in the (x, ξ)-space based on classical orthogonal polynomials of x, like Gegenbauer polynomials, will lead to the same conclusions as presented in this section.This is because any monomial of x can be expressed by a finite number of classical orthogonal polynomials.

Artificial neural network basis
In the case of ANNs, universal approximation theorems (see e.g.Ref. [57]) for arbitrary width or depth of the neural network ensure that a large enough network is able to accurately represent any continuous function on a compact set.The graphical representation of an exemplary network used to construct valid GPD models is shown in Fig. 2.
(1) The exemplary network is made out of four layers: input and output layers ((1) and (4), respectively), and two hidden layers ((2) and ( 3)).The signal is processed through the network as follows: i ) The neuron in the input layer receives the value of x and distributes it to all neurons in the first hidden layer (second layer of the network: (2)) via a set of connections.A single weight is associated to each of those connections that is denoted by w 1,i , where in our example i = 1, . . ., 6. ii ) Each neuron in the hidden and output layers processes the incoming signal according to this equation: where o is the weight associated to the connection between i-th neuron from (l − 1) layer and k-th neuron from (l) layer, and b (l) k is the bias parameter.We note that: iii ) To reproduce the form of Eq. ( 11) we use the linear activation function for the neuron in the output layer: where (•) denotes the function argument, and we introduce the ξ-dependence via weights associated to the connections linking the last hidden layer with the output layer: We also set the bias in the output neuron to: This choice means that f j (x) associated to ξ j , see Eq. ( 11), is effectively described by an ANN with a single hidden layer corresponding to the first hidden layer of the full network, and the output layer corresponding to one neuron in the second hidden layer of the full network.
We are free to choose the activation function associated to the neurons in the hidden layers.For the purpose of further demonstration in the first hidden layer we use either the sigmoid function, or the rectifier function (ReLU), ϕ where Θ(•) is the Heaviside step function.In the second hidden layer the linear activation function will be used, ϕ This choice allows us to analytically evaluate the Mellin coefficients using: for the sigmoid function Eq. ( 27): or, for the ReLU function Eq. ( 28): k w (1) 1,k w 1,k (n + 1) + s w b (2) where: independently of the choice of the activation function in the first hidden layer: i,k Φ (2) where Li i (x) is the polylogarithm and These coefficients depend linearly on b k and w i,k , so the later can be found from a number of Mellin moments by solving the system of equations similar to that of Eq. ( 18).Here, A j = (A 0,j , . . ., A i,j , . . ., A N,j , 0, 0) T is the same as for the polynomial basis, while 1,1+j/2 , . . ., w i + (−1) n w ( The size of the last hidden layer (third layer of the network: (3)) is N/2 + 1, while to ensure both the polynomiality and vanishing at |x| = 1 the size of the fist hidden layer (second layer of the network: (2)) must be at least N + 2.
In our approach the coefficients of the first hidden layer, w 1,i and b (2) i , are not constrained by the Mellin moments.Values of these coefficients can be taken random, or they can be fixed with further criteria, e.g. to obtain the best possible agreement between the input PDF and the resulting GPD in x-space.The sensitivity on the choice of w (2) i coefficients becomes small for large N.

Demonstration of models
We illustrate the modelling presented in this section with the help of the following pion GPD [55,56]: We fix the truncation parameter N to either 4 or 8, and we evaluate the corresponding A n,j coefficients from H π (x, ξ).We use those coefficients to reconstruct H π (x, ξ) using either the polynomial or ANN basis (with either the sigmoid or ReLU activation functions).This way we check how well one may reconstruct the underlying GPD only knowing a number of its Mellin moments.The result is shown in Fig. 3 for ξ = 0 and ξ = 0.5.
In this Article we only show the most elementary use of the modelling based on Eq. ( 11).However several possible extensions or modifications are possible.In particular, one may extend the sum in Eq. ( 14) (for polynomial basis) or the number of neurons in the first hidden layer (for ANN basis) to introduce more free parameters that can be "spent" for a better description of the x-dependence and of higher Mellin coefficients for a given power of ξ.One can also try the modelling with an explicit PDF contribution: which automatically fixes all A n,0 coefficients.If needed, vanishing at |x| = |ξ| can be enforced in this way: which typically causes an oscillatory behaviour along the x-axis.One may also introduce a separate D-term [40]: here denoted by D(x/ξ), which will contribute to A n,n+1 coefficients.This separate contribution allows to reproduce GPDs, which are not smooth functions and usually exhibit a kink at |x| = |ξ|.
In our simple demonstration we do not check if the models fulfil the positivity constraint.With the pion model given by Eq. ( 37) used as a benchmark it is particularly difficult.From Eq. ( 3) we see that if the PDF vanishes for x < 0, the corresponding GPD must vanish as well for x < −|ξ|.We are not able to achieve this behaviour easily via our modelling of Mellin moments.However, we note that with a different benchmark model it should be possible to use a numerical enforcement of positivity.This method will be introduced in the next section of this Article.

Principles of modelling
The advantage of modelling GPDs using double distributions is a natural provision of polynomiality by the Radon transform.In addition, a proper reduction to PDFs can be achieved easily, with a welcome possibility of keeping the singularity at x = 0 for only ξ = 0.The GPD support x ∈ [−1, +1] is also naturally encoded in the double distribution support |α| + |β| ≤ 1.
We are modelling the GPD H(x, ξ) with three contributions, each one being the subject to the Radon transform shown in Eq. ( 10): The first term, (1 − x 2 )F C (β, α) gives the "classical" contribution reproducing both the forward limit and the cross-over line ξ = x [58].This term is crucial for the comparison of GPD models with experimental data available today.The second term, (x 2 − ξ 2 )F S (β, α), gives a "shadow" contribution, which vanishes for both ξ = 0 and ξ = x [24].As we will see, its inclusion is important for a proper estimation of model uncertainties when GPDs are constrained by a sparse set of data.The third term, ξF D (β, α), only contributes to the D-term, which is accessible in analyses of amplitudes for exclusive processes (see e.g.Ref. [37]).The inclusion of this term gives an extra flexibility to the model required to reproduce all x-moments of GPDs.
We note that for singlet GPDs, for which even Mellin moments vanish, the prefactors x 2 and ξ 2 in Eq. ( 41) do not break the polynomiality property.Indeed, for odd Mellin moments of order n, double distributions being the subject to the Radon transform shown in Eq. (10) give contributions up to the A n,n−1 coefficient only -since the A n,n coefficient is zero for parity reasons.Therefore, adding either an x 2 or ξ 2 factor only brings contributions up to A n,n+1 , which is compatible with polynomiality.

F C (β, α) contribution
We factorise F C (β, α) into the forward limit, f (β), which is well-known, at least for the GPD H(x, ξ), and the profile function h C (β, α): The prefactor (1 − x 2 )/(1 − β 2 ) arising from Eqs. ( 41) and ( 42) turns out to be convenient to preserve positivity with an ANN-based modelling of h C (β, α).Since in this study we are only interested in the singlet contribution we fix: The profile function must fulfil the following set of properties to ensure the correct behaviour of the GPD in (x, ξ)-space: i ) even parity w.r.t. the β variable, to keep the whole GPD an odd function of x. ii ) even parity w.r.t. the α variable, to keep the whole GPD an even function of ξ, and therefore to hold time reversal symmetry.iii ) the following normalisation: to ensure the proper reduction to the PDF at ξ = 0. iv ) vanishing at the edge of the support region, to avoid any singular behaviour, except at the x = ξ = 0 point, and to help enforcing the positivity property at x ≈ 1.
We use the following model for h C (β, α) fulfilling all the aforementioned requirements: Here, ANN C (|β|, α) is the output of a single artificial neural network.We have chosen a feed forward artificial neural network with a single hidden layer, whose example is shown in Fig. 4. Our choice simplifies the evaluation of the integral in the denominator of Eq. ( 48).The signal is processed by the network in a way similar to the case described in Sect.3.3.Keeping the nomenclature consistent throughout the text, two neurons in the input layer distribute the |β| and α values: while the neuron in the output layer gives: A number of neurons in the hidden layer process information according to this equation: where the biases, b k , and weights, w 1,k , w 2,k , are free parameters of the network, and where ϕ k (•) is the activation function to be fixed by the network architect.For instance, it can be the sigmoid function from Eq. ( 27) or the ReLU function from Eq. (28).For the output layer we have: where we explicitly used the linear activation function.
There is no bias parameter and the weights, w k,1 are the other free parameters of the network.With the sigmoid function for ϕ (2) k (•), the normalisation factor is: The shadow contribution is modelled in a quite similar way to F C (β, α): Here, f (β) has been already given in Eq. (43).The inclusion of this factor may come as a surprise, as the shadow GPD vanishes in the forward limit, but f (β) helps here to fulfil the positivity constraint.Similarly to h C (β, α), h S (β, α) must be an even function with respect to both β and α variables, and we make it to vanish at |α| + |β| = 1.Its normalisation is however different.We require to ensure the vanishing at ξ = 0. We use the following model for h S (β, α), which fulfils all the aforementioned requirements: Here, N S is a scaling parameter, while ANN S (|β|, α) and ANN S (|β|, α) are neural networks.As one may conclude, effectively, we are constructing the shadow contribution by subtracting two GPDs having the same forward limit.One of these GPDs, related to ANN S (|β|, α), will be the subject of modelling.The other one, related to ANN S (|β|, α), provides an arbitrary reference point.For simplicity we take: that is, one of two networks used in F S (β, α) is the same as that used in F C (β, α).In fits presented in Sect.4.6 ANN S (|β|, α) will have the same architecture as ANN S (|β|, α), but different free parameters (weights and biases).

F D (β, α) contribution
This contribution gives the D-term: where: d i are the coefficients of the expansion on the family of Gegenbauer polynomials and N is a truncation parameter.This definition in (x, ξ)-space explicitly gives: where z = x/ξ.

Remarks
In Eq. ( 52) the variable α is scaled by 1 − |β|, such that α/(1 − |β|) spans over the range of (−1, 1).Such a standardisation of variables is typical for neural networks, as it leads to a faster convergence and it allows to describe equally well the dependencies on all input variables.The factorisation expressed by Eq. ( 42) is arbitrary.In principle, the neural network could absorb some, if not all, information about PDFs.For instance, we could express a double distribution F in the following way: where δ and γ are powers driving the behaviour of PDF at x → 0 and x → 1, respectively, and where the normalisation factor, as used in the denominator of Eq. ( 48), is not needed anymore.This kind of factorisation may be useful when (semi-)inclusive and exclusive data will be simultaneously used to constrain GPDs.

Fits to data
We now discuss technicalities of constraining GPD models defined in (β, α)-space from data.The examples presented in this section demonstrate possible applications of our approach, but also its performance in conditions imposed by currently available data.We start with the presentation of results, and then we elaborate more on aspects of our numerical analysis like minimisation, propagation of uncertainties, treatment of outliers, regularisation and positivity enforcement.
In the first test-case we train our model on N pts = 400 points uniformly covering the domain of −4 < log 10 (x) < log 10 (0.95), −4 < log 10 (ξ) < log 10 (0.95), t = 0 and Q 2 = 4 GeV 2 .That is, in this scenario, x = ξ data are used to constrain our ANN model.The purpose of this test is to check if our approach can be used to reproduce such GPD models as GK.All three contributions from Eq. ( 41) are used in the fit, even if we a priori know that GK does not include any D-term.It is because F C (β, α) and F S (β, α) may contribute to the D-term due to the (1 − x 2 ) and (x 2 − ξ 2 ) factors, respectively.With F D (β, α) included in the fit we are able to compensate for those contributions.We do not enforce the positivity constraint, as the GK model violates the inequality shown in Eq. ( 3).Enforcing this positivity inequality would lead to a discrepancy between our model and GK.The result is shown in Fig. 5 for Mellin moments and in Fig. 6 for the x-dependence.We observe a very good agreement between fitted ANN model and GK, hence we conclude to a successful outcome of this test.
In the second test-case only x = ξ data are used to constrain our GPD model.These are N pts = 200 points uniformly covering the domain of −4 < log 10 (x = ξ) < log 10 (0.95), t = 0 and Q 2 = 4 GeV 2 .The purpose of this test is to demonstrate how our approach can be used to reconstruct GPDs from amplitudes for processes like DVCS and TCS, when described at LO.The positivity inequality is not enforced.The result is shown in Fig. 7. Since we are not interested in the D-term (it is not constrained by x = ξ data, but can be accessed elsewhere, namely in dispersive analyses) the contribution coming from F D (β, α) is removed from the presentation of results.The contribution to the D-term generated by F C (β, α) and F S (β, α) is small with respect to the uncertainties and will be discussed in more detail in the next paragraph.In Fig. 7 one may observe exploding uncertainties, except for the ξ = x line.This behaviour is expected.From the same figure one may judge on the importance of including the shadow contribution, F S (β, α).We see that F C (β, α) is over-constrained by the necessity of reproducing the GPD at both ξ = 0 and ξ = x lines.Since h C (β, α) is normalized, see Eqs. (46) and (48), the neural network lacks the flexibility allowing for a significant contribution to the uncertainties in kinematic domains that are unconstrained by data.
The conditions used in the third test-case are the same as in the second one, except now the positivity inequality is enforced to show its impact on the reduction of uncertainties.The result is shown in Fig. 8.We observe a significant reduction of uncertainties due to positivity.The contributions coming from the F C (β, α) and F S (β, α) terms, together with the D-term induced by those two terms, is shown in Fig. 9.Both contributions are substantial.
In this analysis the minimisation, i.e. the procedure of constraining free parameters, is done with a genetic algorithm [59].This algorithm is iterative.In a single iteration many sets of free parameters, referred to in the literature as "candidates", are simultaneously checked against a "fitness function".After evaluation, the best candidates, i.e. those characterised by the best values of the fitness function, are "crossed over" with a hope of obtaining even better candidates to be used in the next iteration.The cross-over is followed by a "mutation", where a number of free parameters is randomly changed, allowing for a significant reduction of the risk of converging to a local minimum.We note that, since in a given iteration of the genetic algorithm the fitness function is simultaneously evaluated for all candidates, multithreading computing can be employed to improve the performance of the minimisation.
The employed fitness function is the root mean squared relative error [60]: Here, the sum runs over N pts points probing (x, ξ) phase-space, while H 0 (x i , ξ i ) and H(x i , ξ i ) are GPDs given by the reference (here GK) and ANN models, respectively.The RMRSE allows to avoid significant differences of contributions to the fitness function coming from various domains of (x, ξ).
The pseudo-data used in this analysis do not have uncertainties.However, we are still interested to know what are the uncertainties of models in domains unconstrained by data, i.e. what is the effect of the sparsity of data, in particular when only x = ξ data are used in the fit and the neural network was designed to explore the whole x = ξ domain.To estimate this kind of uncertainties we repeat the minimization multiple times, each time starting the genetic algorithms with a random set of initial parameters.We refer to the outcome of such a single minimization as "a replica".A replica is a possible realisation of the GPD model.All replicas reproduce the data used in the minimization, but because of the flexibility of ANNs their behaviour in unconstrained kinematic domains can be considered random.We note that this randomness can be unintentionally suppressed if one uses too small networks or restricts too much the values of weights and biases.
We use the spread of 100 replicas to estimate the model uncertainty in a given kinematic point.Because of the complexity of our fits and unavoidable problems in the minimisation, like falling into local minima, it sometimes happens that a single replica gives values that are very different compared to those given by other replicas.This problem is typical in phenomenological studies (see e.g.Ref. [61]), and requires finding and removing the problematic replicas, or suppressing them in the estimation of uncertainties.In this analysis we employ the following algorithm: for a given population of replicas i ) we randomly choose 1000 points of (x, ξ), ii ) for a given point we evaluate the GPD values that the replicas give in that point, and then we check which of those values are outliers by applying the 3σ-rule [62], iii ) the replicas giving values identified as outliers in more than 10% points are entirely removed from the estimation of uncertainties.
As in other analyses using effectively nonparametric models, a regularisation must be used to avoid biased results caused by over-fitting.Without a regularisation, an ANN tends to describe the training data extremely precisely, resulting in a minimal variance evaluated for these data.It does not mean though that the ANN will describe equally well any other data, i.e. that it will have a predictive power.In general, a bias may appear, because too much attention was paid to describe the training data, and the ANN does not describe well the general trends.Many types of regularisation methods exist, and the selection of a given method typically depends on the problem being under consideration.In this analysis we use the so-called dropout method [63].In this method in a given iteration of the minimisation algorithm (here: the genetic algorithm) a predefined fraction of neurons (here: 10%) is randomly dropped, i.e. some neurons become inactive, they do not process the signal.The output of other neurons is correspondingly scaled to compensate the loss.Effectively, in each iteration a different architecture of ANN is probed, avoiding focusing on details only characterising the training sample.
The genetic algorithm allows for a straightforward implementation of penalisation methods to avoid unwanted results, in particular those not fulfilling binarylike conditions, like inequalities.We use this feature to enforce the positivity condition.To achieve this, for each set of free parameters (i.e. for each candidate) we check by how much we can scale F S (β, α) by changing the N S parameter to saturate the condition given by Eq. ( 3).We check this in 10000 points covering the (x, ξ < x) domain.If we are not able to scale F S (β, α) to respect the constraint in all of those points, the candidate is discarded, i.e. it is dropped from the population of candidates.If we can, we scale F S (β, α) accordingly, i.e. we take the highest allowed value of |N S |.This simple method of enforcing the positivity constraints requires a substantial computing power, but gives satisfactory results.We conclude that up to the numerical noise the constraint is fulfilled, which we demonstrate with Fig. 10 showing the x-dependencies of replicas with respect to the positivity bound for few values of ξ.We note that after the minimisation, F C (β, α) and F S (β, α) may not fulfill alone the positivity constraint given by the chosen PDF, but they do in sum.It indicates that the inclusion of F S (β, α) gives an extra freedom to achieve the positivity.
Fig. 5 Result of constraining the ANN model using 400 points evaluated with the GK model [25][26][27] for sea quarks: x = ξ case, positivity not enforced, all three contributions (F C (β, α), F S (β, α), F D (β, α)) used in the presentation of results.The plot shows the first Mellin coefficient for odd moments (even moments strictly vanish) evaluated for the GK (black dots) and ANN (solid bands) models.The lack of GPD flexible parameterisations fulfilling all required theoretical constraints has been an obstacle slowing down the completion of global fits to experimental data comparable to those currently achieved in the PDF community.Often GPD models are either too rigid to accommodate the measurements, or insufficiently constrained from the theoretical point of view.Since the different theoretical requirements are hardly met in common analytic expressions, neural networks are appealing tools to obtain flexible, yet complex, parameterisations.To the best of our knowledge, our present study is the first providing concrete elements to actually build GPD models based on neural network techniques.
In this Article we discuss different strategies to model GPDs in (x, ξ) and (β, α) spaces.In all cases GPD Mellin moments and polynomiality play a central role.We show in particular that a neural network description of double distributions is an efficient way to implement polynomiality, positivity, discrete symmetries, as well as the reduction to any freely chosen PDF in the forward limit.
An increased control of systematic effects is highly desirable in front of the precision era of GPD physics and extractions.To this aim we also address questions of a practical nature when dealing with many parameters: existence of local minima in optimisation routines, presence of outliers in statistical data analysis, etc.
In view of future GPD fits to experimental DVCS data, we pay a special attention to singlet GPDs, and to contributions with both vanishing forward limit and vanishing DVCS amplitude at leading order in perturbation theory.Our parameterisations may thus be used as tools to study the model dependency affecting extractions of GPDs and derived quantities, like GPD Mellin moments and the total angular momentum of hadrons.Since these parameterisations fulfil theory driven constraints, they may be conveniently used in current and future analyses of GPDs, or in connection e.g. to lattice QCD, either in x-space or through Mellin moments.

Fig. 2
Fig. 2 Scheme of an exemplary artificial neural network used to represent a single GPD.The example is for the truncation parameter N = 4.

Fig. 3
Fig.3The pion model[55,56] (black solid line) given by Eq. (37) reconstructed from either five (N = 4, top row) or nine (N = 8, bottom row) of its first Mellin moments as a function of x for either ξ = 0 (left column) or ξ = 0.5 (right column).The reconstruction is done using a polynomial basis (red dashed line) or an ANN basis with either sigmoid activation function (blue dotted line) or ReLU activation function (green dash-dotted line).

Fig. 4
Fig. 4 Scheme of an exemplary artificial neural network used to represent the profile function in the double distribution model.

Fig. 6 Fig. 7
Fig.6Result of constraining the ANN model using 400 points evaluated with the GK model[25][26][27] for sea quarks: x = ξ case, positivity not enforced, all three contributions (F C (β, α), F S (β, α), F D (β, α)) shown.See the text for more details.The comparison is for (left) ξ = x, (center) ξ = 0.1 and (right) ξ = 0.5.The dashed lines denote the GK model, while the bands represent the result of the fit in the form of a 68% confidence level.

Fig. 8
Fig. 8 Result of constraining the ANN model using 200 points evaluated with the GK model [25-27] for sea quarks: x = ξ case, positivity enforced, F D (β, α) not shown.See the text for more details.The comparison is for (left) x = ξ, (center) ξ = 0.1 and (right) ξ = 0.5.The dashed lines denote the GK model, while the bands represent the result of the fit in the form of a 68% confidence level.The inner bands show F C (β, α) contribution, alone, while the outer bands are for F C (β, α) + F S (β, α).The regions excluded by the positivity constraint are denoted by the hatched bands.