Teaching to extract spectral densities from lattice correlators to a broad audience of learning-machines

We present a new supervised deep-learning approach to the problem of the extraction of smeared spectral densities from Euclidean lattice correlators. A distinctive feature of our method is a model-independent training strategy that we implement by parametrizing the training sets over a functional space spanned by Chebyshev polynomials. The other distinctive feature is a reliable estimate of the systematic uncertainties that we achieve by introducing several ensembles of machines, the broad audience of the title. By training an ensemble of machines with the same number of neurons over training sets of fixed dimensions and complexity, we manage to provide a reliable estimate of the systematic errors by studying numerically the asymptotic limits of infinitely large networks and training sets. The method has been validated on a very large set of random mock data and also in the case of lattice QCD data. We extracted the strange-strange connected contribution to the smeared $R$-ratio from a lattice QCD correlator produced by the ETM Collaboration and compared the results of the new method with the ones previously obtained with the HLT method by finding a remarkably good agreement between the two totally unrelated approaches.


I. INTRODUCTION
The problem of the extraction of hadronic spectral densities from Euclidean correlators, computed from numerical lattice QCD simulations, has attracted a lot of attention since many years (see Refs. , the works on the subject of which we are aware of, and Refs.[31,32] for recent reviews).At zero temperature, the theoretical and phenomenological importance of hadronic spectral densities, strongly emphasized in the context of lattice field theory in Refs.[1,11,13,14,[17][18][19], is associated with the fact that from their knowledge it is possible to extract all the information needed to study the scattering of hadrons and, more generally, their interactions.
From the mathematical perspective, the problem of the extraction of spectral densities from lattice correlators is equivalent to that of an inverse Laplace-transform operation, to be performed numerically by starting from a discrete and finite set of noisy input data.This is a notoriously ill-posed numerical problem that, in the case of lattice field theory correlators, gets even more complicated because lattice simulations have necessarily to be performed on finite volumes where the spectral densities are badly-behaving distributions.
In Ref. [14], together with M. Hansen and A. Lupo, one of the authors of the present paper proposed a method to cope with the problem of the extraction of spectral densities from lattice correlators that allows to take into account the fact that distributions have to be smeared with sufficiently well-behaved test functions.Once smeared, finite volume spectral densities become numerically manageable and the problem of taking their infinite volume 1.By introducing a discrete functional-basis, with elements Bn(E), that is dense in the space of squareintegrable functions f (E) in the interval [E0, ∞) with E0 > 0, any such function can exactly be represented as f (E) = ∞ n=0 cnBn(E).With an infinite number of basis functions (N b = ∞) and by randomly selecting an infinite number (Nρ = ∞) of coefficient vectors c = (c0, • • • , cN b ), one can get any possible spectral density.This is the situation represented by the filled blue disk.If the number of basis functions N b and the number of randomly extracted spectral densities Nρ are both finite one has a training set that is finite and that also depends on N b .This is the situation represented in the first disk on the left.The other two disks schematically represent the situations in which either N b or Nρ is infinite.
limit is mathematically well defined.The method of Ref. [14] (HLT method in short) has been further refined in Ref. [21] where it has been validated by performing very stringent tests within the two-dimensional O(3) non-linear σ-model.
In this paper we present a new method for the extraction of smeared spectral densities from lattice correlators that is based on a supervised deep-learning approach.
The idea of using machine-learning techniques to address the problem of the extraction of spectral densities from lattice correlators is certainly not original (see e.g.Ref. [15,16,[22][23][24][25][26][27][28][29]).The great potential of properlytrained deep neural networks in addressing this problem We generate several training sets Tσ(N b , Nρ) built by considering randomly chosen spectral densities.These are obtained by choosing Nρ random coefficients vectors with N b entries, see FIG. 1.For each spectral density ρ(E) we build the associated correlator C(t) and smeared spectral density ρσ(E), where σ is the smearing parameter.We then distort the correlator C(t), by using the information provided by the statistical variance of the lattice correlator (the one we are going to analyse at the end of the trainings), and obtain the input-output pair (Cnoisy(t), ρσ(E)) that we add to Tσ(N b , Nρ).We then implement different neural networks with Nn neurons and at fixed N = (Nn, N b , Nρ) we introduce an ensemble of machines with Nr replicas.Each machine r = 1, • • • , Nr belonging to the ensemble has the same architecture and, before the training, differs from the other replicas for the initialization parameters.All the replicas are then trained over Tσ(N b , Nρ) and, at the end of the training process, each replica will give a different answer depending upon N .
is pretty evident from the previous works on the subject.These findings strongly motivated us to develop an approach that can be used to obtain trustworthy theoretical predictions.To this end we had to address the following two pivotal questions 1. is it possible to devise a model independent training strategy?
2. if such a strategy is found, is it then possible to quantify reliably, together with the statistical errors, also the unavoidable systematic uncertainties?
The importance of the first question can hardly be underestimated.Under the working assumption, supported by the so-called universal reconstruction theorems (see Refs. [33][34][35]), that a sufficiently large neural network can perform any task, limiting either the size of the network or the information to which it is exposed during the training process means, in fact, limiting its ability to solve the problem in full generality.Addressing the second question makes the difference between providing a possibly efficient but qualitative solution to the problem and providing a scientific numerical tool to be used in order to derive theoretical predictions for phenomenological analyses.
In order to address these two questions, the method that we propose in this paper has been built on two pillars FIG. 3. Flowchart illustrating the three-step procedure that, after the training, we use to extract the final result.Here C(t) represents the input lattice correlator that, coming from a Monte Carlo simulation, is affected by statistical noise.We call Cc(t) the c-th bootstrap sample (or jackknife bin) of the lattice correlator with c = 1, • • • , Nc.In the first step, Cc(t) is fed to all the trained neural networks belonging to the ensemble at fixed N and the corresponding answers ρpred σ (E, N , c, r) are collected.In the second step, by combining in quadrature the widths of the distributions of the answers as a function of the index c (∆ latt σ (E, N )) and of the index r (∆ net σ (E, N )) we estimate the error ∆ stat σ (E, N )) at fixed N .At the end of this step we are left with a collection of results ρpred σ (E, N ) ± ∆ stat σ (E, N ).In the third and last step, the limits N → ∞ are studied numerically and an unbiased estimate of ρpred σ (E) and of its error ∆ tot σ (E), with the latter also taking into account the unavoidable systematics associated with these limits, is finally obtained.
1. the introduction of a functional-basis to parametrize the correlators and the smeared spectral densities of the training sets in a model independent way; 2. the introduction of the ensemble of machines, the broad audience mentioned in the title, to estimate the systematic errors1 .
A bird-eye view of the proposed strategy is given in FIG. 1, FIG. 2 and in FIG. 3. The method is validated by using both mock and true lattice QCD data.In the case of mock data the exact result is known and the validation test is quite stringent.In the case of true lattice QCD data the results obtained with the new method are validated by comparison with the results obtained with the HLT method.The plan of the paper is as follows.In Section II we introduce and discuss the main aspects of the problem.In Section III we illustrate the numerical setup used to obtain the results presented in the following sections.In Section IV we describe the construction of the training sets and the proposed model-independent training strategy.In Section V we illustrate the procedure that we use to extract predictions from our ensembles of trained machines and present the results of a large number of validation tests performed with random mock data.Further validation tests, performed by using mock data generated by starting from physically inspired models, are presented in Section VI.In Section VII we present our results in the case of lattice QCD data and the comparison with the HLT method.We draw our conclusions in Section VIII and discuss some technical details in the two appendixes.

II. THEORETICAL FRAMEWORK
The problem that we want to solve is the extraction of the infinite-volume spectral density ρ(E) from the Euclidean correlators computed numerically in lattice simulations.These are performed by introducing a finite spatial volume L 3 , a finite Euclidean temporal direction T and by discretizing space-time, where a is the so-called lattice spacing and (τ, n) are the integer space-time coordinates in units of a.In order to obtain the infinite-volume spectral density one has to perform different lattice simulations, by progressively increasing L and/or T , and then study numerically the L → ∞ and T → ∞ limits.In the T → ∞ limit the basis functions b T (t, E) become decaying exponentials and the correlator is given by where ρ L (E) = 0 for E < E 0 and the problem is that of performing a numerical inverse Laplace-transform operation, by starting from a discrete and finite set of noisy input data.This is a classical numerical problem, arising in many research fields, and has been thoroughly studied (see e.g.Refs.[37,38] for textbooks discussing the subject from a machine-learning perspective).The problem, as we are now going to discuss, is particularly challenging from the numerical point of view and becomes even more challenging and delicate in our Quantum Field Theory (QFT) context because, according to Wightman's axioms, QFT spectral densities live in the space of tempered distributions and, therefore, cannot in general be considered smooth and well behaved functions.Before discussing though the aspects of the problem that are peculiar to our QFT context, it is instructive to first review the general aspects of the numerical inverse Laplace-transform problem that, in fact, is ill-posed in the sense of Hadamard.To this end, we start by considering the correlator in the infinite L and T limits, and we assume that our knowledge of C(t) is limited to the discrete and finite set of times t = aτ .Moreover we assume that the input data are affected by numerical and/or statistical errors that we call ∆(aτ ).
The main point to be noticed is that, in general, the spectral density ρ(E) contains an infinite amount of information that cannot be extracted from the limited and noisy information contained into the input data C(aτ ).As a consequence, in any numerical approach to the extraction of ρ(E) a discretization of Eq. ( 4) has to be introduced.Once a strategy to discretize the problem has been implemented, the resulting spectral density has then to be interpreted as a "filtered" or (as is more natural to call it in our context) smeared version of the exact spectral density, where K(E, ω) is the so-called smearing kernel, explicitly or implicitly introduced within the chosen numerical approach.
There are two main strategies (and many variants of them) to discretize Eq. ( 4).The one that we will adopt in this paper has been introduced and pioneered by Backus and Gilbert [39] and is based on the introduction of a smearing kernel in the first place.We will call this the "smearing" discretization approach.The other approach, that is more frequently used in the literature and that we will therefore call the "standard" one, is built on the assumption that spectral densities are smooth and well-behaved functions.Before discussing the smearing approach we briefly review the standard one, by putting the emphasis on the fact that also in this case a smearing kernel is implicitly introduced.
In the standard discretization approach the infinitevolume correlator is approximated as a Riemann sum, under the assumption that the infinite-volume spectral density is sufficiently regular to have C(aτ ) − Ĉ(aτ ) ≪ ∆(aτ ) (7) for N E sufficiently large and σ sufficiently small.By introducing the "veilbein matrix" Êτm ≡ e −τ aEm , and the associated "metric matrix" in energy space, Eq. ( 6) is then solved, By using the previous expressions we can now explain why the problem is particularly challenging and in which sense it is numerically ill-posed.On the numerical side, the metric matrix Ĝ is very badly conditioned in the limit of large N E and small σ.Consequently, the coefficients g τ (E n ) become huge in magnitude and oscillating in sign in this limit and even a tiny distortion of the input data gets enormously amplified, In this sense the numerical solution becomes unstable and the problem ill-posed.Another important observation concerning Eqs.(10), usually left implicit, concerns the interpretation of ρ(E n ) as a smeared spectral density.By introducing the smearing kernel and by noticing that, as a matter of fact, the spectral density has to be obtained by using the correlator C(aτ ) (and not its approximation Ĉ(aτ ), to be considered just as a theoretical device introduced in order to formalize the problem), we have Consistency would require that K(E n , ω) = δ(E n − ω) but this cannot happen at finite T and/or N E .In fact K(E n , ω) can be considered a numerical approximation of δ(E n −ω) that, as can easily be understood by noticing that has an intrinsic energy-resolution proportional to the discretization interval σ of the Riemann's sum.A numerical study of K(E n , ω) at fixed E n reveals that for ω ≥ E n the kernel behaves smoothly while it oscillates wildly for ω < E n and small values of σ.A numerical example is provided in FIG. 4.
Once the fact that a smearing operation is unavoidable in any numerical approach to the inverse Laplacetransform problem has been recognized, the smearing discretization approach that we are now going to discuss appears more natural.Indeed, the starting point of the smearing approach is precisely the introduction of a smearing kernel and this allows to cope with the problem also in the case, relevant for our QFT applications, in which spectral densities are distributions.
By reversing the logic of the standard approach, that led us to Eq. ( 12), in the smearing approach the problem is discretized by representing the possible smearing kernels as functionals of the coefficients g τ according to Notice that we are now considering the correlator C LT (aτ ) at finite T and L and this allows to analyse the results of a single lattice simulation.The main observation is that in the T → ∞ limit any target kernel can exactly be represented as a linear combination of the b ∞ (aτ, ω) = exp(−aωτ ) basis functions (that are indeed dense in the functional-space L 2 [E 0 , ∞] of squareintegrable functions).With a finite number of lattice times, the smearing kernel is defined as the best possible approximation of K target (E, ω), i.e. the coefficients g τ (E) are determined by the condition Once the coefficients are given, the smeared spectral density is readily computed by relying on the linearity of the problem, where now In the smearing approach one has Moreover, provided that ρ(E) exists, it can be obtained by choosing as the target kernel a smooth approximation to δ(E − ω) depending upon a resolution parameter σ, e.g.
and by considering the limits lim in the specified order.When σ is very small, the coefficients g τ (E) determined by using Eq. ( 17) become gigantic in magnitude and oscillating in sign, as in the case of the coefficients obtained by using Eq (10).In fact the numerical problem is ill-posed independently of the approach used to discretize it.We now come to the aspects of the problem that are peculiar to our QFT context.Hadronic spectral densities are the Fourier transforms of Wightman's functions in Minkowski space, where P = ( Ĥ, P ) is the QCD four-momentum operator and the Ôi 's are generic hadronic operators.According to Wightman's axioms, vacuum expectation values of field operators are tempered distributions.This implies that also their Fourier transforms, the spectral densities, are tempered distributions in energy-momentum space.It is therefore impossible, in general, to compute numerically a spectral density.On the other hand, it is always possible to compute smeared spectral densities, where the kernels K i (p) are Schwartz functions.In this paper, to simplify the discussion and the notation, we focus on the dependence upon a single spacetime variable, where the operators ÔI and ÔF might also depend upon other coordinates.The spectral density associated with W (x) is given by and, to further simplify the notation, we don't show explicitly the dependence of ρ(E) w.r.t. the fixed spatial momentum p.The very same spectral density appears in the two-point Euclidean correlator for positive Euclidean time t.In the last line of the previous equation we have used the fact that, because of the presence of the energy Dirac-delta function appearing in Eq. ( 27), the spectral density vanishes for E < E 0 where E 0 is the smallest energy that a state propagating between the two hadronic operators ÔI and ÔF can have.On a finite volume the spectrum of the Hamiltonian is discrete and, consequently, the finite-volume spectral density ρ L (E) becomes a particularly wild distribution, the sum of isolated Dirac-delta peaks in correspondence of the eigenvalues E n (L) of the finite volume Hamiltonian.By using the previous expression one has that and this explains why we have discussed the standard approach to the discretization of the inverse Laplacetransform problem by first taking the infinite-volume limit.Indeed if the Riemann's discretization of Eq. ( 6) is applied to C L (aτ ) and all the E m 's are different from all the E n (L)'s one gets ĈL (aτ ) = 0! On the one hand, also in infinite volume, spectral densities have to be handled with the care required to cope with tempered distributions and this excludes the option of using the standard approach to discretize the problem in the first place.On the other hand, in some specific cases it might be conceivable to assume that the infinite volume spectral density is sufficiently smooth to attempt using the standard discretization approach.This requires that the infinite-volume limit of the correlator has been numerically taken and that the associated systematic uncertainty has been properly quantified (a non-trivial numerical task at small Euclidean times where the errors on the lattice correlators are tiny).
The smearing discretization approach can be used either in finite-and infinite-volume and the associated systematic errors can reliably be quantified by studying numerically the limits of Eq. (22).For this reason, as in the case of the HLT method of Ref. [14], the target of the machine-learning method that we are now going to discuss in details is the extraction of smeared spectral densities.

III. NUMERICAL SETUP
In order to implement the strategy sketched in FIG. 1, FIG. 2, FIG. 3 and discussed in details in the following sections, the numerical setup needs to be specified.In this section we describe the data layout that we used to represent the input correlators and the output smeared spectral densities and then provide the details of the architectures of the neural networks that we used in our study.

A. Data layout
In this work we have considered both mock and real lattice QCD data and have chosen the algorithmic parameters in order to be able to extract the hadronic spectral density from the lattice QCD correlator described in Section VII.Since in the case of the lattice QCD correlator the basis function b T (t, E) (see Eq. ( 1)) is given by with T = 64a, we considered the same setup also in the case of mock data.These are built by starting from a model unsmeared spectral density ρ(E) and by computing the associated correlator according to and the associated smeared spectral density according to Since in the case of the lattice QCD spectral density, in order to compare the results obtained with the proposed new method with the ones previously obtained in Ref. [40], we considered a Gaussian smearing kernel, also in the case of mock data we made the choice We stress that there is no algorithmic reason behind the choice of the Gaussian as the smearing kernel and that any other kernel can easily be implemented within the proposed strategy.Among the many computational paradigms available within the machine-learning framework, we opted for the most direct one and represented both the correlator and the smeared spectral density as finite dimensional vectors, that we used respectively as input and output of our neural networks.More precisely, the dimension of the input correlator vector has been fixed to N T = 64, coinciding with the available number of Euclidean lattice times in the case of the lattice QCD correlator.The inputs of the neural networks are thus the 64-components vectors C = {C(a), C(2a), • • • , C(64a)}.The output vectors are instead given by ρσ = {ρ σ (E min ), • • • , ρσ (E max )}.As in Ref. [40], we have chosen to measure energies in units of the muon mass, m µ = 0.10566 GeV, and set E min = m µ and E max = 24m µ .The interval [E min , E max ] has been discretized in steps of size m µ /2.With these choices our output smeared spectral densities are the vectors ρσ with N E = 47 elements corresponding to energies ranging from about 100 MeV to 2.5 GeV.
The noise-to-signal ratio in (generic) lattice QCD correlators increases exponentially at large Euclidean times.For this reason it might be numerically convenient to choose N T < T /a and discard the correlator at large times where the noise-to-signal ratio is bigger than one.According to our experience with the HLT method, that inherits the original Backus-Gilbert regularization mechanism, it is numerically inconvenient to discard part of the information available on the correlator provided that the information contained in the noise is used to conveniently regularize the numerical problem.Also in this new method, as we are going to explain in subsection IV B, we use the information contained in the noise of the lattice correlator during the training process and this, as shown in Appendix B, allows us to conveniently use all the available information on the lattice correlator, i.e. to set N T = T /a, in order to extract the smeared spectral density.
We treat σ, the width of the smearing Gaussian, as a fixed parameter by including in the corresponding training sets only spectral functions that are smeared with the chosen value of σ.This is a rather demanding numerical strategy because in order to change σ the neural networks have to be trained anew, by replacing the smeared spectral densities in the training sets with those corresponding to the new value of σ.Architectures that give the possibility to take into account a variable input parameter, and a corresponding variable output at fixed input vector, have been extensively studied in the machine learning literature and we leave a numerical investigation of this option to future work on the subject.In this work we considered two different values, σ = 0.44 GeV and σ = 0.63 GeV, that correspond respectively to the smallest and largest values used in Ref. [40].

B. Architectures
By reading the discussion on the data layout presented in the previous subsection from the machine-learning point of view, we are in fact implementing neural networks to solve a R N T → R N E regression problem with N T = 64 and N E = 47 which, from now on, fix the dimension of the input and output layers of the neural networks.
There are no general rules to prefer a given network architecture among the different possibilities that have been considered within the machine-learning literature It consist of three 1D convolutional layers with an increasing number of maps followed by two fully connected layers.The two blocks are intermediated by one flatten layer.The column denoted by "Size" reports the shape of the signal produced by the corresponding layer.The stride of the filters is set to 2 in such a way that the dimension of the signal is halved at each 1D convolutional layer thus favouring the neural network to learn a more abstract, and possibly more effective, representation of the input data.As activation functions we use the LeakyReLu with negative slope coefficient set to −0.2.The neurons with activation functions are also provided with biases.The output is devoid of activation function in order not to limit the output range.The bottom line reports the total number of trainable parameters.
and it is common practice to make the choice by taking into account the details of the problem at hand.For our analysis, after having performed a comparative study at (almost) fixed number of parameters of the so-called Multilayer perceptron and convolutional neural networks, we used feed-forward convolutional neural networks based on the LeNet architecture introduced in Ref. [41].
We studied in details the dependence of the output of the neural networks upon their size N n and, to this end, we implemented three architectures that we called arcS, arcM and arcL.These architectures, described in full details in TABLEs I, II and III, differ only for the number of maps in the convolutional layers.The number of maps are chosen so that the number of parameters of arcS:arcM:arcL are approximately in the proportion 1 : 2 : 3.For the implementation and the training we employed both Keras [42] and TensorFlow [43].

IV. MODEL INDEPENDENT TRAINING
In the supervised deep-learning framework a neural network is trained over a training set which is representative of the problem that has to be solved.In our case the inputs to each neural network are the correlators C and the target outputs are the associated smeared spectral densities ρσ .As discussed in the Introduction, our main goal is to devise a model-independent training strategy.To this end, the challenge is that of building a training set which contains enough variability so that, once trained, the network is able to provide the correct answer, within the quoted errors, for any possible input correlator.
As a matter of fact, the situation in which the neural network can exactly reconstruct any possible function is merely ideal.That would be possible only in absence of errors on the input data and with a neural network with an infinite number of neurons, trained on an infinitely large and complex training set.This is obviously impossible and in fact our goal is the realistic task of getting an output for the smeared spectral density as close as possible to the exact one by trading the unavoidable limited abilities of the neural network with a reliable estimate of the systematic error.In order to face this challenge we used the algorithmic strategy described in FIG. 1, FIG. 2 and in FIG. 3. In our strategy, • the fact that the network cannot be infinitely large is parametrized by the fact that N n (the number of neurons) is finite; • the fact that during the training a network cannot be exposed to any possible spectral density is parametrized by the fact that N b (the number of basis functions) and N ρ (the number of spectral densities to which a network is exposed during the training) are finite (see FIG. 1); • the fact that at fixed the answer of a network cannot be exact, and therefore has to be associated with an error, is taken into account by introducing an ensemble of machines, with N r replicas, and by estimating this error by studying the distribution of the different N r answers in the N r → ∞ limit (see FIG. 2); • once the network (and statistical) errors at fixed N are given, we can study numerically the N → ∞ limits and also quantify, reliably, the additional systematic errors associated with these unavoidable extrapolations (see FIG. 3).
We are now going to provide the details concerning the choice of the functional basis that we used to parametrize the spectral densities and to build our training sets.

A. The functional-basis
In our strategy we envisage studying numerically the limit N b → ∞ and, therefore, provided that the systematic errors associated with this extrapolation are properly taken into account, there is no reason to prefer a particular basis w.r.t.any other.For our numerical study we used the Chebyshev polynomials of the first kind as basis functions (see for example Ref. [44]).
The Chebyshev polynomials T n (x) are defined for x ∈ [−1, 1] and satisfy the orthogonality relations In order to use them as a basis for the spectral densities that live in the energy domain E ∈ [E 0 , ∞), we introduced the exponential map and set Notice that the subtraction of the constant term T n (x(E 0 )) has been introduced in order to be able to cope with the fact that hadronic spectral densities vanish below a threshold energy E 0 ≥ 0 that we consider an unknown of the problem.With this choice, the unsmeared spectral densities that we use to build our training sets are written as and vanish identically for E ≤ E 0 .Once E 0 and the coefficients c n that define ρ(E; N b ) are given, the correlator and the smeared spectral density associated with ρ(E; N b ) can be calculated by using Eq. ( 32) and Eq. ( 33) 2 .Each ρ(E; N b ) that we used in order to populate our training sets has been obtained by choosing E 0 randomly in the interval [0.2, 1.3] GeV with a uniform distribution and by inserting in Eq. ( 39) the coefficients where the r n 's are N b uniformly distributed random numbers in the interval [−1, 1] and ε is a non-negative parameter that we set to 10 −7 .Notice that with this choice of the coefficients c n the Chebyshev series of Eq. ( 39) is convergent in the N b → ∞ limit and that the resulting spectral densities can be negative and in general change sign several times in the interval E ∈ [0, ∞).Notice also that up to a normalization constant, that it will turn to be irrelevant in our training strategy in which the input data are standardized as explained in subsection IV C, the choice of the interval [−1, 1] for the r n random numbers is general.
A few examples of unsmeared spectral densities generated according to Eq. 39 are shown in FIG. 5.As it can be seen, with the choice of the coefficients c n of Eq. ( 40), the larger is the number of terms in the Chebyshev series of Eq. ( 39) the richer is the behaviour of the resulting unsmeared spectral densities in terms of local structures.This is an important observation and a desired feature.Indeed, a natural concern about the choice of Chebyshev 2 By representing also the smearing kernels and the basis functions on a Chebyshev basis, the orthogonality relations of Eq. ( 36) can conveniently be exploited to speedup this step of the numerical calculations.polynomials as basis functions is the fact that we are thus sampling the space of regular functions while, as we pointed out in section II, on a finite volume the lattice QCD spectral density is expected to be a discrete sum of Dirac-delta peaks (see Eq. ( 29)).Here we are relying on the fact that the inputs of our networks will be Euclidean correlators, that in fact are smeared spectral densities, and will also be asked to provide as output the smeared spectral densities ρσ (E).This allows us to assume that, provided that the energy smearing parameter σ and N b are sufficiently large, the networks will be exposed to sufficiently general training sets to be able to extract the correct smeared spectral densities within the quoted errors.To illustrate this point we consider in FIG.6 a continuous spectral density ρ(E) and approximate its smeared version through a Riemann's sum according to where The previous few lines of algebra show that the smearing of a continuous function ρ(ω) can be written as the smearing of the distribution defined in Eq. ( 42), a prototype of the finite-volume distributions of Eq. ( 29), plus the approximation error Σ σ,∆E (E) = ρσ (E) − ρδ,σ (E).
The quantity Σ σ,∆E (E), depending upon the spacing ∆E of the Dirac-delta peaks and the smearing parameter σ, is expected to be sizeable when σ ≤ ∆E and to become irrelevant in the limit σ ≫ ∆E.A quantitative example is provided in FIG.6 where Σ σ,∆E (E) is shown at fixed ∆E for two values of σ.In light of this observation, corroborated by the extensive numerical analysis that we have performed at the end of the training sessions to validate our method (see subsection V B and, in particular, FIG.16), we consider justified the choice of Chebyshev polynomials as basis functions provided that, as is the case in this work, the energy smearing parameter σ is chosen sufficiently large to not be able to resolve the discrete structure of the finite-volume spectrum.

B. Building the training sets
Having provided all the details concerning the functional basis that we use to parametrize the unsmeared spectral densities, we can now explain in details the procedure that we used to build our training sets.
A training set T σ (N b , N ρ ) contains N ρ input-output pairs.Each pair is obtained by starting from a random unsmeared spectral density, parametrized at fixed N b according to Eq. ( 39) and generated as explained in the previous subsection.Given the unsmeared spectral density ρ(E; N b ), we then compute the corresponding correlator vector C (by using Eq. ( 32)) and, for the two values of σ that we used in this study, the smeared spectral densities vectors ρσ (by using Eq. ( 33)).From each pair (C, ρσ ) we then generate an element (C noisy , ρσ ) of the training set at fixed N b and σ, that we obtain, as we are now going to explain, by distorting the correlator C using the information provided by the noise of the lattice correlator C latt (see Section VII).
In order to cope with the presence of noise in the input data that have to be processed at the end of the training, it is extremely useful (if not necessary) to teach to the networks during the training to distinguish the physical content of the input data from noisy fluctuations.This is particularly important when dealing with lattice QCD correlators for which, as discussed in subsection III A, the noise-to-signal ratio grows exponentially for increasing Euclidean times (see FIG. 7).A strategy to cope with this problem, commonly employed in the neural network literature, is to add Gaussian noise to the input data used in the training.There are several examples in the literature where neural networks are shown to be able to learn rather efficiently by employing this strategy of data corruption (see the already cited textbooks Refs.[37,38]).According to our experience, it is crucially important that the structure of the noise used to distort the training input data resembles as much as possible that of the true input data.In fact, it is rather difficult to model the noise structure generated in Monte Carlo simulations and, in particular, the off-diagonal elements of the covariance matrices of lattice correlators.In the light of these observations we decided to use the covariance matrix Σlatt of the lattice correlator C latt that we are going to analyse in Section VII to obtain C noisy from C.More precisely, given a correlator C, we generate C noisy by extracting a random correlator vector from the multivariate Gaussian distribution C. Data pre-processing A major impact in the machine-learning performance is played by the way the data are presented to the neural network.The learning turns to be more difficult when the input variables have different scales and this is exactly the case of Euclidean lattice correlators since the components of the input vectors decrease exponentially fast.The risk in this case is that the components small in magnitude are given less importance during the training or that some of the weights of the network become very large in magnitude, thus generating instabilities in the training process.We found that a considerable improvement in the validation loss minimization is obtained by implementing the standardization of the input correlators (see next subsection and the central panel of FIG.8) and therefore used this procedure.
The standardization procedure of the input consists in rescaling the data in such a way that all the components of the input correlator vectors in the training set are distributed with average 0 and variance 1.For a given training set T σ (N b , N ρ ) we calculate the N T -component vectors µ and γ whose components are and Each correlator in the training is then replaced by where C noisy is the distorted version of C discussed in the previous subsection.Notice that the vectors µ and γ are determined from the training set before including the noise.Although the pre-processed correlators look quite different from the original ones, since the components are no longer exponential decaying, the statistical correlation in time is preserved by the standardization procedure.At the end of the training, the correlators fed into the neural network for validation or prediction have also to be pre-processed by using the same vectors µ and γ used in the training.

D. Training an ensemble of machines
Given a machine with N n neurons we train it over the training set T σ (N b , N ρ ) by using as loss function the Mean Absolute Error (MAE) where ρpred,i σ (w) is the output of the neural network in correspondence of the input correlator C i noisy .In Eq. ( 48) we used the norm | ρj | and we have explicitly shown the dependence of the predicted spectral density ρpred,i σ (w) upon the weights w of the network.At the beginning of each training session each weight w n (with n = 1, • • • , N n ) is extracted from a Gaussian distribution with zero mean value and variance 0.05.To end the training procedure we rely on the early stopping criterion: the training set is split into two subsets containing respectively the 80% and 20% of the entries of the training set T σ (N b , N ρ ).The larger subset is used to update the weights of the neural network with the gradient descent algorithm.The smaller subset is the socalled validation set and we use it to monitor the training process.At the end of each epoch we calculate the loss function of Eq. ( 48) for the validation set, the so-called validation loss, and stop the training when the drop in the validation loss is less than 10 −5 for 15 consecutive epochs.In the trainings performed in our analysis this occurs typically between epoch 150 and 200.The early stopping criterion provides an automatic way to prevent the neural network from overfitting the input data.
We implement a Mini-Batch Gradient Descent algorithm, with Batch Size (BS) set to 32, by using the Adam optimizer [45] combined with a learning rate decaying according to where e ∈ N is the epoch and η(0) = 2 × 10 −4 is the starting value.The step-function is included so that the learning rate is unchanged during the first 25 epochs.Although a learning rate scheduler is not strictly mandatory, since the Adam optimizer already includes adaptive learning rates, we found that it provides an improvement in the convergence with less noisy fluctuations in the validation loss (see the top panel of FIG. 8).Concerning the BS, we tested the neural network performance by starting from BS=512 and by halving it up to BS=16 (see bottom panel of FIG. 9).Although the performance improves as BS decreases we set BS= 32 in order to cope with the unavoidable slowing down of the training for smaller values of BS.
As we already stressed several times, at fixed N the answer of a neural network cannot be exact.In order to be able to study numerically the N → ∞ limits, the error associated with the limited abilities of the networks  densities, and therefore from the same pairs (C, ρσ ) (see Eq. ( 43)), but with different noisy input correlator vectors C noisy .
In FIG. 9 we show the validation loss as a function of N ρ for the different values of N b , the three different architectures and the two values of σ considered in this study.Each point in the figure has been obtained by studying the distribution of the validation loss at the end of the training within the corresponding ensemble of machines and by using respectively the mean and the standard deviation of each distribution as central value and error.The figure provides numerical evidences concerning the facts that • networks with a finite number N n of neurons, initialized with different weights and exposed to training sets containing a finite amount of information, provide different answers and this is a source of er-rors that have to be quantified; • the validation loss decreases for larger N n because larger networks are able to assimilate a larger amount of information; • the validation loss decreases for larger N ρ since, if N n is sufficiently large, the networks learn more from larger and more general training sets; • at fixed N n and N ρ the networks perform better at smaller values of N b because the training sets exhibit less complex features and learning them is easier.
In order to populate the plots in FIG. 9 we considered several values of N ρ and N b and this is rather demanding from the computational point of view.The training sets that we found to be strictly necessary to quote our final results are listed in TABLE IV.

V. ESTIMATING THE TOTAL ERROR AND VALIDATION TESTS
Having explained in the previous section the procedure that we used to train our ensembles of machines, we can now explain in details the procedure that we use to obtain our results.We start by discussing the procedure that we use to estimate the total error, taking into account both statistical and systematics uncertainties, and then illustrate the validation tests that we have performed in order to assess the reliability of this procedure.

A. Procedure to estimate the total error
The procedure that we use to quote our results for smeared spectral densities by using the trained ensembles of machines, illustrated in FIG. 3 • by computing the standard deviation over the ensemble of machines of ρpred σ (E, N , r), we obtain an estimate of the error, that we call ∆ net σ (E, N ), associated with the fact that at fixed N the answer of a network cannot be exact; • both ∆ latt σ (E, N ) and ∆ net σ (E, N ) have a statistical origin.The former, ∆ latt σ (E, N ), comes from the limited statistics of the lattice Monte Carlo simulation while the latter, ∆ net σ (E, N ), comes from the statistical procedure that we used to populate our training sets T σ (N b , N ρ ) and to train our ensembles of machines at fixed N b .We sum them in quadrature and obtain an estimate of the statistical error at fixed N , • we then study numerically the N → ∞ limits by using a data-driven procedure.We quote as central value and statistical error of our final result where N max , given in Table IV, is the vector with the largest components among the vectors N considered in this study; • in order to check the numerical convergence of the N → ∞ limits to estimate the associated systematic uncertainties ∆ X σ (E), where we define the reference vectors N ref X listed in Table IV.We then define the pull variables and then we weight the absolute value of the difference ρpred with the Gaussian probability that this is not due to a statistical fluctuation, (56) • the total error that we associate to our final result ρpred σ (E) is finally obtained according to (57)   56) are maximised on our setup and, consequently, Eq. ( 57) provides a conservative estimate of the total error.Some remarks are in order here.We don't have a theoretical understanding neither of the dependence of the statistical errors ∆ latt σ (E, N ) and ∆ net σ (E, N ) upon E and N (see Appendix C for a numerical investigation) nor of the rates of convergence of the N → ∞ limits.Gaining this theoretical understanding is a task that goes far beyond the scope of this paper.The procedure that we have devised to quote our results, that at first sight might appear too complicated, has therefore to be viewed as just one of the possible ways to perform conservative plateaux-analyses of the N → ∞ limits.This explains our choice of the points N ref X given in Table IV that, our numerical setup, provide the most conservative estimates of the systematic errors.The stringent validation tests that we have performed with mock data, and that we are going to discuss below, provide a quantitative evidence that the results obtained by implementing this procedure can in fact be trustworthy used in phenomenological applications.

B. Validation
In order to illustrate how the method described in the previous subsection can be applied in practice, we now consider an unsmeared spectral density ρ(E) that has not been used in any of the trainings and that has been obtained as described in subsection IV A, i.e. by starting from Eq. ( 39), but this time with N b = 1024, i.e. a number of basis functions which is twice as large as the largest N b employed in the training sessions.
From this ρ(E) we then calculate the associated correlator C and the smeared spectral density corresponding to σ = 0.44 GeV.We refer to the true smeared spectral density as ρtrue σ (E), while we call ρpred σ (E) the final predicted result.Both the unsmeared spectral density ρ(E) (grey dotted curve, partially visible) and the expected exact result ρtrue σ (E) (solid black curve) are shown in the top panel of FIG.12.
Starting from the exact correlator C corresponding to ρ(E) we have then generated N c = 800 bootstrap samples from the distribution of Eq. ( 44), thus simulating the outcome of a lattice Monte Carlo simulation.The N c samples have been fed into each trained neural network and we collected the answers ρpred σ (E, N , c, r) that, as   shown in FIG. 10 for a selection of the considered values of E, we have then analysed to obtain ∆ latt σ (E, N ) and ∆ net σ (E, N ).The next step is now the numerical study of the limits N → ∞ that we illustrate in FIG.11 for the same values of E considered in FIG.10.The limit N ρ → ∞ (left Relative error budget (%) one in each plot.The red points correspond to the results ρpred that we use to estimate the systematic uncertainties ∆ X σ (E).The blue band correspond to our estimate of the total error ∆ tot σ (E).In the top panel of FIG. 12 we show the comparison of our final results ρpred σ (E) ± ∆ tot σ (E) (blue points) with the true smeared spectral density ρtrue σ (E) (black curve) that in this case is exactly known.The central panel in FIG. 12 shows the relative error budget as a function of the energy.As it can be seen, the systematics errors represent a sizeable and important fraction of the total error, particularly at large energies.The bottom panel of FIG. 12 shows the pull variable and, as it can be seen, by using the proposed procedure to estimate the final results and their errors, no significant deviations from the true result have been observed in this particular case.In order to validate our method we repeated the test just described 2000 times.We have generated random spectral densities, not used in the trainings, by starting again from Eq. ( 39) and by selecting random values for E 0 and N b respectively in the intervals [0.3, 1.3] GeV and [8,1024].
The values of p σ (E) for all the energies and for all the samples is collected in one set, p σ for short, whose normalized distribution is shown in the left panels of FIG. 13 (top panel for σ = 0.44 GeV, bottom panel for σ = 0.63 GeV).The distributions are nicely bell-shaped around 0 and the right panels in FIG. 13 show the comparison between the normalized cumulative distribution functions of p σ with the normal distribution obtained from the mean and variance of p σ .The p-value calculated from the Kolmogorov-Smirnov test is much less than 0.05 in both the cases and therefore the distribution of p σ cannot be considered as a normal one.On the other hand, given the procedure that we use to quote the total error, we observe deviations ρpred σ (E) − ρtrue σ (E) smaller than 2 standard deviations in 95% of the cases and smaller than 3 standard deviations in 99% of the cases.The fact that deviations smaller than 1 standard deviation occur in ∼80% of the cases for σ = 0.44 GeV and ∼85% of the cases for σ = 0.63 GeV (to be compared with the expected 68% in the case of a normal distribution) is an indication of the fact that our estimate of the total error is indeed conservative.Moreover, from FIG. 13 it is also evident that our ensembles of neural networks are able to generalize very efficiently outside the training set.
In order to perform another stringent validation test of the proposed method, we also considered unsmeared spectral densities that cannot be written on the Chebyshev basis.We repeated the analysis described in the previous paragraph for a set of spectral densities generated according to where 0 By plugging Eq. (59) into Eq.( 32) and Eq. ( 33) we have calculated the correlator and smeared spectral density associated to each unsmeared ρ(E).We have generated 2000 unsmeared spectral densities with N peaks = 5000.The position E n of each peak has been set by drawing independent random numbers uniformly distributed in the interval [E 0 , 15 GeV] while E 0 has been randomly chosen in the interval [0.3, 1.3] GeV.The coefficients c n have also been generated randomly in the interval [−0.01, 0.01].
The 2000 trains of isolated Dirac-delta peaks are repre- sentative of unsmeared spectral densities that might arise in the study of finite volume lattice correlators.These are rather wild and irregular objects that our neural networks have not seen during the trainings.FIG 14 shows the plots equivalent to those of FIG 13 for this new validation set.As it can be seen, the distributions of p σ (E) are basically unchanged, an additional reassuring evidence of the ability of our ensembles of neural networks to generalize very efficiently outside the training set and, more importantly, on the robustness of the procedure that we use to estimate the errors.

VI. RESULTS FOR MOCK DATA INSPIRED BY PHYSICAL MODELS
In the light of the results of the previous section, providing a solid quantitative evidence of the robustness and reliability of the proposed method, we investigate in this section the performances of our trained ensembles of machines in the case of mock spectral densities coming from physical models.A few more validation tests, unrelated to physical models, are discussed in Appendix A.
For each test discussed in the following subsections we define a model spectral density ρ(E) and then use Eq.(32) and Eq. ( 33) to calculate the associated exact correlator and true smeared spectral density.We then generate N c = 800 bootstrap samples C c (t), by sampling the distribution of Eq. ( 44), and quote our final results by using the procedure described in details in the previous section.

A. A resonance and a multi-particle plateaux
The first class of physically motivated models that we investigate is that of unsmeared spectral densities exhibiting the structures corresponding to an isolated resonance and a multi-particle plateau.To build mock spectral densities belonging to this class we used the same Gaussian mixture model used in Ref. [24] and given by The sigmoid function θ(E, δ 1 , ζ 1 ) with δ 1 = 0.1 GeV and ζ 1 = 0.01 GeV dumps ρ GM (E) at low energies while the other sigmoid θ(E, δ 2 , ζ 2 ), with ζ 2 = 0.1 GeV, connects the resonance to the continuum part whose threshold is δ 2 .The parameter C 0 = 1 regulates the height of the continuum plateaux which also coincides with the asymptotic behaviour of the spectral density, i.e. ρ GM (E) → C 0 for E → ∞.We have generated three different spectral densities with this model that, in the following, we call ρ GM 1 (E), ρ GM 2 (E) and ρ GM 3 (E) and whose parameters are given in TABLE V.The predicted smeared spectral densities for smearing widths σ = 0.44 GeV and σ = 0.63 GeV are compared to the true ones in FIG. 15.In all cases the predicted result agrees with the true one, within the quoted errors, for all the explored values of E. The quality of the reconstruction of the smeared spectral densities is excellent for E < 1.5 GeV while, at higher energies, the quoted error is sufficiently large to account for the deviation of ρpred σ (E) from ρtrue σ (E).This is particularly evident in the case of ρ GM 1 (E) shown in the top panel.
The same class of physical models can also be investigated by starting from unsmeared spectral densities that might arise in finite volume calculations by considering The parameter E peak parametrizes the position of an isolated peak while the multi-particle part is introduced with the other peaks located at E peak < E 1 ≤ • • • ≤ E N peaks .We have generated two unsmeared spectral densities, ρ peak 1 (E) and ρ peak 2 (E) by using this model.In both cases we set N peaks = 10000, selected random values for the E n 's up to 50 GeV and random values for the c n 's in the range [0, 0.01].In the case of ρ peak 1 (E) we set E peak = 0.8 GeV, C peak = 1 and E 1 = 1.2 GeV.In the case of ρ peak 2 (E) we set instead E peak = 0.7 GeV, C peak = 1 and E 1 = 1.5 GeV.
The predicted smeared spectral densities are compared with the true ones in FIG.16.In both cases ρpred σ (E) is in excellent agreement with ρtrue σ (E).It is worth emphasizing once again that the model in Eq. ( 62) cannot be represented by using the Chebyshev basis of Eq. ( 39) and, therefore, it is totally unrelated to the smooth unsmeared spectral densities that we used to populate the training sets.

B. O(3) non-linear σ-model
In this subsection we consider a model for the unsmeared spectral density that has already been investigated by using the HLT method in Ref. [21].More precisely, we consider the two-particles contribution to the vector-vector spectral density in the the 1+1 dimensional O(3) non-linear σ-model (see Ref. [21] for more details) given by where E th is the two-particle threshold.We considered three mock unsmeared spectral densities, that we call ρ (E) and ρ O(3) 3 (E), differing for the position of the multi-particle threshold, which has been set respectively to E th = 0.5 GeV, 1 GeV and 1.5 GeV.The reconstructed smeared spectral densities for σ = 0.44 GeV and σ = 0.63 GeV are compared with the exact ones in FIG. 17.
The predicted smeared spectral densities are in remarkably good agreement with the true ones in the full energy range and in all cases.This result can be read as an indication of the fact that the smoothness of the underlying unsmeared spectral density plays a crucial rôle in the precision that one can get on ρpred σ (E), also if the problem is approached by using a neural network approach.Indeed, this fact had already been observed and exploited in Ref. [21], where the authors managed to perform the σ → 0 limit of ρO (3) σ (E) with controlled systematic errors by using the HLT method.

C. The R-ratio with mock data
The last case that we consider is that of a model spectral density coming from a parametrization of the experimental data of the R-ratio.The R-ratio, denoted by R(E), is defined as the ratio between the inclusive cross section of e + e − → hadrons and e + e − → µ + µ − and plays a crucial rôle in particle physics phenomenology (see e.g.Refs.[40,47]).In the next section we will present results for a contribution to the smeared R-ratio R σ (E) that we obtained by feeding into our ensembles of trained neural networks a lattice QCD correlator that has been already used in Ref. [40] to calculate the same quantity with the HLT method.
Before doing that, however, we wanted also to perform a test with mock data generated by starting from the parametrization of R(E) given in Ref. [46].This parametrization, that we use as model unsmeared spectral density by setting ρ(E) ≡ R(E), is meant to reproduce the experimental measurements of R(E) for energies E < 1.1 GeV, i.e. in the low-energy region where there  In this section we use our trained ensembles of neural networks to extract the smeared spectral density from a real lattice QCD correlator.
We have considered a lattice correlator, measured by the ETMC on the ensemble described in TABLE VI, that has already been used in Ref. [40] to extract the so-called strange-strange connected contribution to the smeared R-ratio by using the HLT method of Ref. [14].The choice is motivated by the fact that in this case the exact answer for the smeared spectral density is not known and we are going to compare our results Rpred σ (E) with the ones, that we call RHLT σ (E), obtained in Ref. [40].In QCD the R-ratio can be extracted from the lattice FIG.19.The results obtained with our ensembles of neural networks in the case of a real lattice QCD correlator are compared with those obtained with the HLT method (grey, Ref. [14,40]).The top-panel shows the two determinations of the strange-strange connected contribution to the smeared R-ratio for σ = 0.44 GeV while the case σ = 0.63 GeV is shown in the bottom-panel. correlator where J µ (x) is the hadronic electromagnetic current.The previous formula (in which we have neglected the term proportional to e −(T −t)ω that vanishes in the T → ∞ limit) explains the choice that we made in Eq. ( 32) to represent all the correlators that we used as inputs to our neural networks.In Ref. [40] all the connected and disconnected contributions to C latt (t), coming from the fact that J µ = f q f ψf γ µ ψ f with f = {u, d, s, c, b, t}, q u,c,t = 2/3 and q d,s,b = −1/3, have been taken into account.Here we consider only the connected strangestrange contribution, i.e. we set J µ = − 1 3 sγ µ s and neglect the contribution associated with the fermionicdisconnected Wick contraction.From C latt (t) we have calculated the covariance matrix Σlatt that we used to inject noise in all the training sets that we have built

and used in this work (see Section IV B).
Although, as already stressed, in this case we don't know the exact smeared spectral density, we do expect from phenomenology to see a sizeable contribution to the unsmeared spectral density coming from the ϕ(1020) resonance.Therefore, the shape of the smeared spectral density is expected to be similar to those shown in FIG.16 for which the results obtained from our ensembles of neural networks turned out to be reliable.The comparison of Rpred σ (E) and RHLT σ (E) is shown in FIG.19 and, as it can be seen, the agreement between the two determinations is remarkably good.Some remarks are in order here.The HLT method and the new method that we propose here are radically different and the fact that the results obtained with the two procedures are in such a good agreement, especially for E < 1.5 GeV where both methods provide very small errors, is at the same time reassuring and encouraging in view of future applications of spectral reconstruction techniques.
Concerning the errors, there is no significant evidence that one method has better performances than the other.This can be better appreciated in FIG.20 where the relative error is shown in both cases as a function of the energy.
We want to stress that in the new proposed method, as in the HLT case, no prior information on the expected spectral density has to be provided and, therefore, the results for Rpred σ (E) shown in FIG.19 have to be considered as first-principles model-independent determinations of the smeared strange-strange contribution to the R-ratio that we managed to obtain by using supervised deep-learning techniques.
We leave to future work on the subject the task of analysing all the contributions to R(E), as done in Ref. [40], in order to obtain a machine-learning determination of Rσ (E) to be compared with experiments.

VIII. CONCLUSIONS
In this work we have proposed a new method to extract smeared hadronic spectral densities from lattice correlators.The method has been built by using supervised deep-learning techniques and is characterized by the distinctive features to implement a model-independent training strategy and to provide a reliable estimate of both the statistical and systematic uncertainties.
We managed to implement a model-independent training strategy by introducing a basis in the functional space from which we extract the spectral densities that we use to populate our training sets.In order to obtain a reliable estimate of the systematic errors, we introduced the ensembles of machines.
We have shown that, by studying the distribution of the answers of the different machines belonging to an ensemble, at fixed and finite number of neurons, dimension and complexity of the training set, it is possible to quantify reliably the systematic errors associated with the proposed method.To do that, we presented a large number of stringent validation tests, performed with mock data, providing quantitative evidence of the reliability of the procedure that we use to estimate the total error on our final results (see FIG. 13 and FIG.14).
In addition to mock data, we have applied the new proposed method also in the case of a lattice QCD correlator obtained from a simulation performed by the ETM Collaboration.We have extracted the strange-strange connected contribution to the smeared R-ratio and compared the predictions obtained by using our ensembles of trained machines with the ones previously obtained by using the HLT method [14,40].We found a remarkably good agreement between the results obtained by using the two totally unrelated methods that provide total errors of the same order of magnitude.
The proposed method requires the training of many machines with different architectures and dimensions of the training sets.Admittedly, although this is a task that can be completed in a few days on a modern desktop computer, the procedure might end up to be computationally demanding.On the one hand, we don't exclude the possibility that the proposed method can be simplified in order to speed up the required computations.At the same time, on the basis of our experience, we are firmly convinced that a careful study of the different sources of systematic uncertainties is mandatory when dealing with machine-learning techniques and when the aim is to compare theoretical predictions with experiments.In fact, the computational cost of the proposed method is the price that we had to pay for the reliability of the results.
As a final remark we want to stress that in this paper we have taught a lesson to a broad audience of learningmachines.The subject of the lesson, the extraction of smeared spectral densities from lattice correlators, is just a particular topic.The idea of teaching systematically to a broad audience of machines is much more general and can be used to estimate reliably the systematic errors in many other applications.errors that are much larger in these cases for σ = 0.44 GeV.This had to be expected and is a desired feature.Indeed, the inverse problem becomes harder when one tries to reconstruct sharp peaks with high resolution.The presence of statistical noise in the input data prevents any method to provide a very precise result, especially at high energies.Our neural network strategy, in these exceptional situations, provides a large error and this is a further reassuring evidence of the robustness of the proposed procedure.In Section III we emphasized that the noise-to-signal ratio of generic lattice correlators grows exponentially with the Euclidean time (see FIG. 7) and that for this reason it might be convenient to reduce the number N T of components of the input vector C.
The approaches based on the Backus-Gilbert regularization method have a built-in mechanism to suppress noisy input data and by using these methods one can safely feed the whole information contained in the correlator into the algorithm.In general, this cannot be expected to be true when supervised deep-learning approaches are employed.Even though neural networks have been proven to be robust in handling noisy data, too much noise can have a bad impact on the performances since a big fraction of the effort in the minimization algorithm is put in the suppression of the noise and in learning the distinction between outliers and effective information rather than in learning general features of the problem.
In the following we investigate the dependence of the training performances, obtained by injecting the noise of the lattice correlator in the training sets as explained in subsection IV B, upon N T .In FIG.22 we show the validation loss as a function of N T .As it can be seen, the performances of the networks improve as N T increases but the validation loss reaches a saturation point around N sat T = 40.On the one hand, N sat T can be considered as the maximum number of time slices of the correlator from which meaningful physical information can be extracted.On the other hand, despite the corruption of the input data for N T > N sat T , including more time slices does not compromise the quality of the trainings and we could safely set N T = T /a = 64.In illustrating our method we have stressed that there are two sources of statistical errors that contribute to ∆ stat σ (E, N ) defined in Eq. ( 52): ∆ latt σ (E, N ), coming from the Monte Carlo simulation, and ∆ net σ (E, N ), coming from our ensemble of trained machines.In this appendix we study the dependence of ∆ latt  3).Therefore, even small statistical fluctuations in the correlator turn into enhanced errors in the spectral density at high energy.At fixed energy the Monte Carlo error ∆ latt σ does not show a significant dependence upon N .This is a reassuring evidence that ∆ latt σ and ∆ net σ are two decorrelated sources of uncertainties.Conversely, at fixed energy ∆ net σ decreases by increasing N ρ and/or N n , thus confirming our working hypothesis that an infinitely large neural network, trained over an infinitely large dataset, is able to provide the exact answer independently of the initialization parameters.∆ net σ is instead stable with respect to changes in the dimension N b of the functional-basis.This stability can be ascribed to the fact that we are targeting the extraction of smeared spectral densities at sufficiently large values of σ so that the local complexity introduced in the unsmeared spectral densities by increasing N b is washed out.A trend in ∆ net σ as a function of N b may thus become visible when considering smaller values of σ.

4 FIG. 4 .
FIG.4.Top panel : Smearing kernel of Eq. 12 at En = 5 for three values of σ.The reconstruction is performed by setting E0 = 1 and Emax = 10 as UV cutoff in Eq. 5 and by working in lattice units with a = 1.The kernel is smooth for ω ≥ En while it presents huge oscillations for ω < En.In the case σ = 0.1 these oscillations range from −10 126 to +10 121 .Bottom panel : Corresponding coefficients gτ (En) (see Eq. 10) in absolute value for the first 20 discrete times.Extendedprecision arithmetic is mandatory to invert the matrix Ĝnm.This test, in which det( Ĝ) = O 10 −15700 for σ = 0.1, has been performed by using 600-digits arithmetic

FIG. 5 .
FIG.5.Examples of unsmeared spectral densities generated according to Eq. 39 using the Chebyshev polynomial as basis functions for different N b .For all the examples we set E0 = 0.8 GeV.As it can be seen, by generating the cn coefficients according to Eq. (40), the larger the number of terms of the Chebyshev series the richer is the behaviour of the spectral density in terms of local structures.
FIG.6.The black dashed curve is a generic continuous unsmeared spectral density ρ(E) while the vertical dashed lines indicate the energy sampling defining the support of ρ δ (ω) according to Eq. (42).The blue and red curves are the smeared version of respectively ρ(ω) and ρ δ (ω).The green curve represents Σσ,∆E(E) = ρσ(E) − ρδ,σ (E).The spacing is set to ∆E = 0.1 GeV.The smearing kernel is a normalized Gaussian function of width σ and as resolution width we refer to the full width at half maximum 2 √ 2 ln 2 σ .Top-panel : the resolution width is equal to ∆E and it is such that the single peaks are resolved.Bottom-panel : in this case 2 √ 2 ln 2 σ ≫ ∆E and the two smeared functions are almost undistinguishable.

2 .
C as mean vector and as covariance the matrix Σlatt normalized by the factor C(a) C latt (a) In order to be able to perform a numerical study of the limits N → ∞, we have generated with the procedure just described several training sets, corresponding to N b = {16, 32, 128, 512} and N ρ = {50, 100, 200, 400, 800} × 10 3 .At fixed N b , each training set includes the smaller ones, i.e. the training set T σ (N b , 100 × 10 3 ) includes the training set T σ (N b , 50 × 10 3 ) enlarged with 50 × 10 3 new samples.

FIG. 8 .
FIG. 8. Numerical tests to optimize the training performances.The validation loss functions of the trainings of the Nr = 20 replica machines belonging to the same ensemble are plotted with the same color.All the trainings refer to arcL, Nρ = 50 × 10 3 , N b = 512 and σ = 0.44 GeV.If not otherwise specified the data are pre-processed, BS=32 and the learning rate scheduler is implemented.Top panel : Comparison of the validation losses with and without the learning scheduler rate defined in Eq. (49).Central panel : Comparison of the validation losses by implementing or not the input data pre-processing procedure described in subsection IV C. Bottom panel : Comparison of the validations losses at different values of the BS parameter.

3 × 10 − 3 4 × 10 − 3 6 × 10 FIG. 9 .
FIG. 9. Validation losses at the end of the trainings as functions of the training set dimension Nρ.The top-plot corresponds to σ = 0.44 GeV and the bottom-plot to σ = 0.63 GeV.The different colors correspond to the different architectures and different markers to the different number of basis functions N b considered in this study.The error on each point correspond to the standard deviation of the distribution of the validation loss in an ensemble of Nr = 20 machines.As expected, the validation loss decreases by increasing Nρ (more general training) and/or Nn (larger and therefore smarter networks).Moreover, at fixed Nn and Nρ the neural network performs better at smaller values of N b as consequence of the fact that the training set exhibits less complex features and learning them is easier.
, is the following • given an ensemble of N r machines trained at fixed N , we feed in each machine r belonging to the ensemble all the different bootstrap samples (or jackknife bins) C c (aτ ) of the input correlator (c = 1, • • • , N c ) and obtain a collection of results ρpred σ (E, N , c, r) for the smeared spectral density that depend upon N , c and r; • at fixed N and r we compute the bootstrap (or jackknife) central values ρpred σ (E, N , r) and statistical errors ∆ latt σ (E, N , r) and average them over the ensemble of machines, ρpred σ (E, N ) = 1 N r Nr r=1 ρpred σ (E, N , r) ,

2 5 10 FIG. 10 .
FIG. 10.Each row corresponds to a different energy, as written on the top of each plot, and in all cases we set σ = 0.44 GeV and N = N max = (arcL, 512, 800 × 10 3 ).Left panels: the x-axes corresponds to the replica index r = 1, • • • , Nr = 20 running over the entries of the ensemble of machines.The plotted points correspond to the results ρpred σ (E, N , r) ± ∆ latt σ (E, N , r) obtained by computing the bootstrap averages and errors of the c = 1, • • • , Nc = 800 results ρpred σ (E, N , c, r).Right panels: distributions of the central values ρpred σ (E, N , r).The means of these distributions correspond to the results ρpred σ (E, N ) that are shown in the left panels as solid blue lines.The widths correspond to the network errors ∆ net σ (E, N ) that are represented in the left panels with two dashed blue lines.The blue bands in the left panels correspond to ∆ stat σ (E, N ), defined in Eq. (52).

FIG. 11 .
FIG. 11.Numerical study of the N → limits.Each row corresponds to a different energy, as written on the top of each plot, and all cases we set σ = 0.44 GeV.Left panels: Study of the limit Nρ → ∞ at fixed N b = 512 and Nn =arcL.Central panels: Study of the limit Nn → ∞ at fixed N b = 512 and Nρ = 800 × 10 3 .Right panels: Study of the limit N b → ∞ at fixed Nn =arcL and Nρ = 800 × 10 3 .All panels: The horizontal blue line is the final central value, corresponding to N = N max = (arcL, 512, 800 × 10 3 ) (blue points) and it is common to all the panels within the same row.The blue band represents instead ∆ tot σ (E) calculated by the combination in quadrature of all the errors (see Eq. 57).The red points correspond to N ref X and the grey ones to the other trainings, see TABLE IV.

FIG. 12 .
FIG.12.Final results for σ = 0.44 GeV.Top panel : The final results predicted by the neural networks (blue points) is compared with the true result (solid black).The dotted grey curve in the background represents the unsmeared spectral density.Central panel : Relative error budget as a function of the energy.Bottom panel : Deviation of the predicted result from the true one in units of the standard deviation, see Eq. (58).

FIG. 13 .
FIG.13.Left panel : Normalized distribution of the significance defined in Eq. (58).The distribution is calculated over 2000 validation samples generated on the Chebyshev functional space.Right panel : comparison of the normalized cumulative distribution functions (CDF) of the observed distribution of pσ (blue) and of the normal one obtained from the mean and variance of pσ (red).The black arrow represents the Kolmogorov-Smirnov statistics (DKS), i.e. the magnitude and the position of the maximum deviation between the two CDFs in absolute value.

FIG. 14 .
FIG.14.The as FIG.13but for 2000 samples of unsmeared spectral densities generated according to Eq. (59) and thus totally unrelated to the data used during the training.

FIG. 16 .
FIG.16.Reconstructed smeared spectral densities (points with errors) compared to the true ones (solid lines) for two unsmeared spectral densities simulating a finite volume distribution, see Eq. (62).The unsmeared ρ(E) is represented by vertical lines located in correspondence of E = En with heights proportional to cn.

FIG. 17 .
FIG. 17. Predicted smeared spectral densities for the twoparticles contribution to the vector-vector correlator in the 1 + 1 dimensional O(3) non-linear σ model for different values of the two-particle threshold: 0.5 GeV (top panel), 1 GeV (central panel) and 1.5 GeV (third panel).The solid lines and the dashed curve correspond respectively to ρpred σ (E) and ρ O(3) (E).

FIG. 20 .
FIG.20.Relative error on the smeared spectral density obtained from our ensembles of neural networks (colored points) and the HLT method (grey points) for σ = 0.44 GeV (top panel) and σ = 0.63 GeV (bottom panel).

FIG. 21 .
FIG.21.Predicted smeared spectral densities for the three models of unsmeared spectral densities given in Eq. (A1).The solid lines correspond to ρtrue σ (E).The vertical grey lines are in correspondence of the isolated Dirac-delta peaks and the height is proportional to the associated weights.

FIG. 22 .
FIG. 22. Validation losses at the end of the training for different number of time slices in the input data and for different training set size Nρ.This test is performed with N b = 512 and σ = 0.44 GeV.Each point with the associated error bar comes from the average of Nr = 20 trainings (see Section IV).

FIG. 23 .
FIG. 23.Behaviour of ∆ latt σ (E, N ) (blue) and ∆ net σ (E, N ) (red) with respect to N = {Nρ, Nn, N b } at four different values of the energy.The data refer to σ = 0.44 GeV and have been obtained by analysing the 2000 samples built on the Chebyshev basis already used in Section V B. The central values and the associated error bars come from the average over the 2000 samples.

σ
and ∆ net σ upon the energy E and upon N .The investigation is carried out by collecting ∆ latt σ and ∆ net σ for the 2000 mock samples built on the Chebyshev basis and used in Section V B to validate our procedure.For each combination of N the 2000 results are averaged and the standard deviation is calculated.The results are shown in FIG.23 at σ = 0.44 GeV and at four different energies spanning the output energy range.We observe that both ∆ latt σ and ∆ net σ increase in magnitude as the energy increases.This phenomenology is a well-known feature of the inverse Laplace transform problem in which the contributions to the correlation function coming from high energies are strongly suppressed by the exponential basis (see Eq. (

TABLE I .
arcS : the smallest neural network architecture used in this work.The architecture is of the type feed-forward and the structure can be read from top to bottom of the table.

TABLE II .
arcM : the medium-size architecture used in this work.See TABLE I for the description.

TABLE III .
arcL: the largest architecture used in this work.See TABLE I for the description.
FIG.7.Noise-to-signal ratio of the lattice correlator C latt (t), discussed in Section VII, whose covariance matrix is used to inject the statistical noise in the training sets.

TABLE IV .
List of the vectors N used to estimate the systematic errors associated with the N → ∞ limits.With these choices the errors ∆ X σ (E) defined in Eq. (

TABLE V .
Additional parameters used to generate the mock unsmeared spectral densities based on a Gaussian mixture model, Eq. 60.The corresponding functions are plotted in FIG. 15 (dashed grey).
are two dominant structures associated with the mixed ρ and ω resonances, a rather broad peak at E ≃ 0.7 GeV, and the narrow resonance ϕ(1020).Given this rich structure, shown as the dashed grey curve in FIG.18, this is a much more challenging validation test w.r.t. the one of the O(3) model discussed in the previous subsection.σ (E) − Rtrue σ (E) doesn't exceed the quoted error for all the considered values of E. VII.THE R-RATIO WITH LATTICE QCD DATA