Interpolating amplitudes

Motivation

• Numerical multi-loop amplitudes are valuable beyond 2 → 2 scattering, but their lengthy evaluation at millions of phase-space points makes direct Monte Carlo use unfeasible.Analytic calculations become rapidly impractical as external legs, scales, or masses increase, while numerical methods scale more tractably with the number of scales.Source: §1 · p. 3

Approach

• Interpolation of precomputed amplitude values is presented as a promising first step toward making numerically known, non-analytic multi-loop amplitudes accessible for phenomenology.The approach evaluates amplitudes at selected points and interpolates them at phase-space points generated by Monte Carlo programs.Source: §1 · p. 3

Scope

• Whereas existing grid-based applications mainly address 2 → 2 amplitudes with two kinematic variables, 2 → 3 interpolation requires accurate approximation over a five-dimensional phase space.The higher-dimensional setting makes achieving precision highly non-trivial.Source: §1 · p. 3

Contribution

• The paper studies state-of-the-art multidimensional interpolation methods to connect mathematical interpolation research with physics use cases that have largely remained below five dimensions.Its stated purpose is to assess practical methods for loop-amplitude functions whose specific properties may not be captured by standard general function spaces.Source: §1 · p. 3

2.1 What is an amplitude

• Virtual amplitudes become especially costly when additional kinematic scales or masses prevent fast Monte Carlo evaluation of analytic corrections, motivating loop-amplitude interpolation.Tree-level and one-loop amplitudes are used as test cases while targeting two-loop applications.Source: §2.1 · pp. 4–5
• Phase spaces are parameterized by d variables in the unit hypercube, with five dimensions for 2 →3 processes and two dimensions for 2 →2 processes.Common parameter choices include energies, Mandelstam variables, and angular variables.Source: §2.1 · p. 5

2.2 Our goal

• The interpolation goal is to approximate an expensive amplitude component from as few evaluations as possible while remaining sufficiently close to the original function.The expensive component is modeled as a function from the d-dimensional unit hypercube to the real numbers.Source: §2.2 · p. 5

2.3 How to define the approximation error

• Approximation error is defined using the L1 distance between differential cross sections, making it variable-independent and weighting regions by their contribution to the total cross section.This choice guarantees bounded regional cross-section errors, but prioritizes total-cross-section precision over relative precision in low-contribution tails.Source: §2.3 · p. 5

2.4 How to prioritize different phase-space regions

• The weight w can guide interpolation by transforming the target, concentrating regression points, refining adaptive regions, or applying a variable transformation; these approaches are studied for each method.The paper does not discuss these weighting strategies further as standalone techniques, and limited testing found no major gains from the transformation-based approaches.Source: §2.4 · p. 6

2.5 The test functions

• The study spans five test functions, mainly five-dimensional amplitudes for q̄q →t̄tH and gg →t̄tH, plus a two-dimensional gg →Hg amplitude.The test set includes integrable Coulomb-type singularities that are subtracted and a divergent gg →Hg region excluded through a transverse-momentum cutoff of approximately 25 GeV.Source: §2.5 · pp. 6, 7, 9, 10

2.6 How to use test function symmetries

• Interpolation methods exploit discrete test-function symmetries through construction constraints, duplicated data, or reduced domains, and every method uses at least one such technique.These choices reduce redundant evaluations or restrict interpolation to a symmetry-reduced region.Source: §2.6 · pp. 10–11

2.7 How to evaluate the approximation error

• Low-discrepancy Sobol sampling provides only marginally better error convergence than random Monte Carlo sampling, with errors stabilizing after about 1000 testing samples.This comparison uses sparse-grid interpolations of f2 with varying data-point counts.Source: §2.7 · p. 11, Figure 1 · p. 12

3 Polynomial interpolation

• The method constructs a degree-(n−1) polynomial passing through exactly n interpolation nodes in the univariate case.Source: §3 · p. 13
• The interpolant is evaluated in barycentric form because this purely polynomial representation is numerically stable and simple.The barycentric formula can also express rational interpolation through suitable weights, but the weights used here correspond to polynomial interpolation.Source: §3 · p. 13
• Polynomial approximation error depends on the node locations, motivating node choices that minimize the relevant product over the domain.This choice avoids requiring precise prior knowledge of the nth derivative f^(n).Source: §3 · p. 13
• Equidistant nodes cause the Runge phenomenon, with approximation errors worsening near boundaries as polynomial degree increases.The passage illustrates this behavior for 10 nodes on [−1;1].Source: §3 · p. 13

Chebyshev nodes of the first kind

• Chebyshev nodes of the first kind minimize the node product in the infinity-norm sense on [−1;1].Source: Chebyshev nodes of the first kind · p. 13
• These nodes achieve a uniform result for the relevant interpolation-error control over the interval.Source: Chebyshev nodes of the first kind · p. 13
• Although nonlinear rational methods such as AAA [44] are highly efficient in one dimension, their multidimensional generalizations are not developed well enough for this study.Source: Chebyshev nodes of the first kind · p. 13

Chebyshev polynomials

• Chebyshev polynomials of the first kind have the interpolation nodes as zeros and are orthogonal under a specified weight.Source: Chebyshev polynomials · p. 14
• Chebyshev nodes are equidistant in φ = arccos(x), making the polynomial expansion an even Fourier series in φ.Source: Chebyshev polynomials · p. 14
• The transformation from function values to Chebyshev coefficients is a discrete cosine transform.Source: Chebyshev polynomials · p. 14
• For n ≳40, rounding errors make the monomial and direct coefficient forms unsuitable, so the barycentric interpolation form is preferred.Source: Chebyshev polynomials · p. 14

Chebyshev nodes of the second kind

• Chebyshev nodes of the second kind avoid the Runge phenomenon by placing points at extrema and endpoints rather than at zeros of T_n(x).Source: Chebyshev nodes of the second kind · p. 14
• Their nested structure lets adaptive interpolation reuse results from smaller grids contained in larger grids, although point spacing is nonuniform.Source: Chebyshev nodes of the second kind · p. 14
• Endpoint inclusion can be impractical for scattering amplitudes because evaluation at x_2 = 0 may lose precision and evaluation near x_1 = 1 is typically slower.Source: Chebyshev nodes of the second kind · p. 14

Gauss nodes

• Gauss nodes are defined by Legendre-polynomial zeros, while Gauss–Lobatto nodes are defined by their extrema.Both are closely related to Chebyshev nodes and are often considered superior.Source: Gauss nodes · p. 14
• Gauss quadrature is exact for polynomials through degree 2n −1, compared with degree n −1 for Clenshaw–Curtis quadrature on Chebyshev nodes.Despite this formal advantage, the practical approximation errors are very close.Source: Gauss nodes · p. 15
• Because Gauss nodes are substantially harder to compute and provide similar practical errors, the paper does not consider them further.Source: Gauss nodes · p. 15

3.1 Approximation error scaling

• At Chebyshev nodes, interpolation is logarithmically close to the best polynomial interpolant of the same degree.The approximation error depends on the smoothness of f [42, 48, 49].Source: §3.1 · p. 15

3.2 Multiple dimensions

• For multiple dimensions, the construction uses an outer tensor product of one-dimensional Chebyshev nodes, allowing different node counts n_i in each dimension.The resulting interpolant can be formed by nested one-dimensional interpolation or a tensor-product decomposition.Source: §3.2 · p. 15
• Tensor-product interpolation converges progressively more slowly as dimension d increases, an effect known as the curse of dimensionality.The multidimensional error scaling is qualitatively analogous to the one-dimensional estimates, with n replaced by the per-dimension counts n_i, whose relevant scale grows with d.Source: §3.2 · p. 15

3.3 Dimensionally adaptive grid

• Dimensionally adaptive sampling can tune per-dimension node counts, with the reported ratio n1:n2:n3:n4:n5 ≈ 2:4:2:3:1.The fifth dimension has faster coefficient decay, while discrete symmetries in the third, fourth, and fifth dimensions create checkerboard patterns of zero coefficients.Source: §3.3 · p. 15
• The best dimension ratio depends on the total number of samples and can only be selected a posteriori, rather than fixed universally.The fixed ratio is therefore a practical approximation to the per-n optimum.Source: §3.3 · p. 17
• Figure 4 compares optimized and non-adaptive polynomial interpolations across sampling schemes, symmetry handling, weighting, and amplitude representations.The “a” method interpolates the amplitude directly, whereas “f” uses Section 2.4 option 1; “half-domain” denotes Section 2.6 option 3, and the remaining methods use option 2.Source: §3.3 · p. 17

3.4 Beyond the full grid

• Even after optimizing cardinal-direction bounds, tensor-product grids can oversample diagonals because diagonal directions are longer than edges.The upper-right regions of the coefficient slices contain significantly smaller coefficients, motivating removal of those regions from the basis.Source: Figure 2 · p. 16, §3.4 · p. 17
• Beyond full grids, the proposed strategy selects a restricted multi-index set and corresponding non-tensor-product nodes, rather than retaining the tensor-product construction.This general framework includes hyperbolic cross, lower, sparse, and adaptively selected polynomial index sets, alongside several specialized node constructions [53,54,55,56,57,58,59,60,61,62,63,64,65].Source: §3.4 · p. 17
• The investigated instance uses an Lp truncation, retaining multi-indices satisfying ||i⃗||p ≤ k.This provides a concrete representative of the broader non-tensor-product interpolation family described in and implemented in.Source: §3.4 · p. 17
• Its nodes are selected greedily in Leja order from the tensor-product node set, with each next point chosen to maximize the stated product criterion.The resulting interpolation performance is evaluated for different choices of p in Figure 5.Source: §3.4 · p. 17

3.5 Effects of noisy data

• Multiplicative noise is used to assess how replacing high-precision evaluations of a with low-precision values affects approximation error.Higher demanded precision makes evaluating a(x⃗) more time-consuming, motivating a practical precision target that ideally does not increase approximation error.Source: §3.5 · p. 17
• Figure 6 reports the results for different σ values, using intentionally high noise levels to make the effects more visible.Source: §3.5 · p. 17

• Polynomial interpolation is a safe and dependable default, with exponential convergence for analytic one-dimensional functions and accurate barycentric evaluation.Source: §3.6 · pp. 17, 20
• In multiple dimensions, tensor-product grids suffer from the curse of dimensionality, while nearby singularities and nonanalytic boundaries worsen convergence.Source: §3.6 · p. 20
• Weighting helps the peaky f5 but slightly hinders other functions; factoring out f5's 1/sin2 θH peak yields a modest but notable improvement.Source: §3.6 · p. 20
• Domain reduction significantly improves f3 and f4, barely affects f1 and f2, and slightly worsens f5, where data duplication performs better.Source: §3.6 · p. 20
• Taking a leading-order ratio does not improve f2 but notably improves f4 at lower precisions; dimensional adaptivity notably helps only f1.Source: §3.6 · p. 20
• Noise sets a lower error boundary, requiring σ ≤ϵtarget, while σ ≲1/3ϵtarget permits the target accuracy with negligible noise effects.Source: §3.6 · p. 20

4.1 B-spline basis functions

• B-spline basis functions are recursively defined on knot vectors, with boundary knot placement ensuring exact polynomial representation throughout the interpolation domain.Uniform knot vectors provide the straightforward construction, while not-a-knot vectors omit selected data-point knots and are better defined for odd degrees.Source: §4.1 · p. 21

4.2 B-spline interpolants

• Multidimensional B-spline interpolants use tensor products and require (n+1)^d interpolation nodes satisfying Schoenberg–Whitney conditions, while retaining numerical stability.The basis functions are non-negative, form a partition of unity, have local support, and can be evaluated efficiently, especially with uniform knots.Source: §4.1 · p. 21, §4.2 · p. 22

4.3 Discussion

• Interpolating the amplitude is better for f1 and f3, whereas including the weight helps f5 with little training data but hurts at high training-data levels.The comparison considers interpolation of the amplitude versus the test function.Source: §4.3 · p. 22
• Tensor-product B-splines inherit polynomial-like flexibility problems, making required sampling hard to predict and potentially excessive for a target precision.Their practical error estimates are therefore optimistic, motivating adaptive sparse-grid extensions.Source: §4.3 · p. 22
• Boundary evaluations are required, which can be problematic when amplitude-evaluation times diverge near boundaries or when shifting the domain leaves phase-space regions uncovered.Extrapolating adaptive sparse-grid bases can alleviate this boundary limitation.Source: §4.3 · p. 22
• B-spline evaluation takes small constant time rather than time proportional to data points, enabling fast surrogates for more precise approximations.A B-spline can approximate an already constructed approximation using cheaply obtained evaluations of that approximation.Source: §4.3 · p. 22
• Quadratic or cubic B-splines significantly outperform linear splines across all test functions; cubic is slightly best for f1 and f2, while quadratic is best for f3, f4, and f5.Figure 8 compares approximation errors for linear, quadratic, and cubic B-splines as training data increase.Source: §4.3 · p. 22, Figure 8 · p. 23

5 Sparse grids

• Classical sparse grids use hierarchical hat-function bases and omit full-grid points that contribute little, thereby alleviating the curse of dimensionality.The section uses hierarchical decomposition because it straightforwardly accommodates spatial adaptivity.Source: §5 · p. 24, §5.1 · p. 24

5.1 Classical sparse grids

• Sparse-grid accuracy decreases only by a logarithmic factor, reaching O(h_k^2(log h_k^-1)^(d−1)), while the construction avoids basis functions highly refined in all directions simultaneously.Source: §5.1 · p. 25

5.2 Boundary treatment

• Boundary support is handled either by adding boundary points, which can erase high-dimensional savings, or by modified basis functions that linearly extrapolate toward boundaries.Sparse grids are best suited when boundary behavior matters less than bulk structure.Source: §5.2 · pp. 25–26

5.3 Spatially adaptive sparse grids

• Spatial adaptivity refines points contributing most to the interpolant, using either largest hierarchical surplus or a support-volume-weighted balanced criterion.The balanced criterion weights surplus by 2^-|l⃗|_1, reducing the risk that refinement becomes trapped in a small region.Source: §5.3 · p. 27

5.5 Discussion

• Balanced refinement is more reliable than greedy refinement, especially for complicated or noisy targets, although greedy refinement can be slightly better for simple functions.Greedy refinement may become trapped in local structures, whereas balanced refinement suppresses local deviations and is therefore preferred for these studies.Source: §5.3 · p. 27, §5.5 · p. 29
• Increasing degree significantly improves all test functions with piecewise polynomial and fundamental B-spline bases, while extended not-a-knot B-splines improve mainly f5.For f3 and f4, cubic and quintic extended not-a-knot splines converge very poorly despite their prior use on high-dimensional tests.Source: §5.4 · p. 27, §5.5 · p. 29
• Ratio interpolation for loop amplitudes significantly improves the low-data regime, but scales worse than direct amplitude interpolation, with higher-degree bases more robust to that scaling.Weighting refinement by the phase-space weight produces only minor differences from balanced refinement for f1 and f3.Source: §5.5 · p. 29
• Spatially adaptive sparse grids provide a significant high-data performance boost through better scaling, despite only modest improvement over non-adaptive full grids at low data.They also support iterative point addition, validation during construction, and stopping at a required point count, unlike non-adaptive methods.Source: §5.5 · pp. 29–30

6 Machine Learning Techniques

• The networks incorporate Lorentz symmetry either through Lorentz-invariant MLP inputs or L-GATr operations constrained to be equivariant under Lorentz transformations.L-GATr uses a geometric-algebra transformer with equivariant linear, attention, normalization, and activation layers, plus a geometric-product MLPBlock.Source: §6 · pp. 33–34, §6.2 · pp. 36–37

6.1 MLP

• MLP(x) uses five phase-space inputs for f1–f4 and two for f5, whereas MLP(s) uses pairwise Lorentz invariants, totaling ten and six inputs respectively.Both versions have five hidden layers with 512 channels; inputs are logarithmically preprocessed and standardized.Source: §6.1 · pp. 34–35

6.2 L-GATr

• L-GATr represents each particle as a 16-component multivector and uses permutation labels to enforce particle-exchange symmetries alongside Lorentz equivariance.Its architecture contains eight attention blocks, 64 multivector channels, 32 scalar channels, and eight attention heads.Source: §6.2 · p. 37

6.3 Discussion

• MLP(s) and L-GATr perform similarly across test functions, while MLP(x) is consistently worse; MLP(s) is slightly better for f5 and low-precision f2, but worse for high-precision f4.These comparisons are shown for training and testing on unweighted samples.Source: §6.3 · p. 37, Figure 16 · p. 39
• Neural-network improvement slows as training data increases, with deterministic data and no large-dataset overfitting indicating training time as the likely bottleneck.Longer MLP training substantially improves performance in the large-data limit, although the limitation merits further study.Source: §6.3 · pp. 37–38
• Machine-learning methods surpass other interpolation approaches at low-to-medium precision for f1–f4, especially the higher-order corrections f2 and f4, but are overtaken at high precision.Neural networks underperform for f5, the only two-dimensional test function studied.Source: §6.3 · p. 38, Figure 16 · p. 39
• Training on uniform or unweighted samples yields similar performance across training-set sizes, allowing naively sampled data to target physically motivated evaluation distributions.Validation samples are unweighted, and the validation set is 10% of training data up to a maximum of 4000 points.Source: §6 · p. 34, §6.3 · p. 38, Figure 17 · p. 40
• Dividing amplitudes by leading-order amplitudes produces a small, stable improvement for f2 and f4, unlike the corresponding polynomial and sparse-grid behavior.Adaptive data addition focused on poorly estimated regions is proposed for future work, potentially combined with uncertainty-estimating ensembles or Bayesian networks.Source: §6.3 · p. 38, Figure 18 · p. 40

Scope and evaluation

• The study benchmarks polynomial, B-spline, spatially adaptive sparse-grid, multilayer-perceptron, and Lorentz-Equivariant Geometric Algebra Transformer interpolation on 5-dimensional t¯tH and 2-dimensional Higgs+jet test functions.The main performance measure is approximation-error scaling with the number of data points, using physically motivated metrics weighted by phase-space contributions to the total cross section.Source: §7 · p. 41

Overall findings

• Machine learning significantly outperforms classical methods for 5-dimensional test functions at 10−2 to 10−4 approximation error, while sparse grids are best for very high precision in lower dimensions.Machine-learning performance stagnates beyond a precision threshold, whereas classical methods continue to scale smoothly; the threshold remains high enough for machine learning to be the practical high-dimensional choice.Source: §7 · p. 41

Training targets

• Training on the ratio a/aLO improves interpolation of higher-order corrections, especially for neural networks and, for some methods, in the low-data regime.Neural-network approaches benefit from ratio training at all training sizes.Source: §7 · p. 41

Limitations

• Two-loop examples were excluded because useful 2 → 3 cases would require training-data generation beyond the paper’s scope.Simple two-loop cases have low-dimensional phase spaces where any interpolation method would perform well.Source: §7 · p. 41

Extendibility

• Adaptive methods and neural networks support adding training data after validation, with sparse grids offering the most granular error-guided refinement.Sparse grids can add single points and identify poorly approximated regions, while neural networks tolerate varied point distributions but require retraining.Source: §7 · p. 41

Resources

• Implementations of the interpolation methods and benchmark test functions are provided in a public GitHub repository.The repository is https://github.com/OlssonA/interpolating_amplitudes.Source: §7 · p. 41

A.1 Phase space parametrisation for pp →t ¯t H

• The pp →t ¯t H test functions f1 to f4 use a 2 →3 phase-space parametrisation with three angular variables and two energy variables.This parametrisation is designed for producing a top-quark pair and a Higgs boson.Source: §A.1 · p. 42
• The angular variables include the Higgs polar angle θH relative to the beam axis and the top-pair solid-angle element dΩt¯t = sinθt dθt dϕt.Source: §A.1 · p. 42
• The Källén function λ(a, b, c) = a2 + b2 + c2 −2ab −2bc −2ca enters the phase-space variables.Its arguments depend on the kinematic variables introduced in the parametrisation.Source: §A.1 · p. 42
• The parametrisation uses the partonic centre-of-mass invariant ŝ = (p3 + p4 + p5)2 and the top-pair invariant st¯t = (p3 + p4)2.Source: §A.1 · p. 42
• The t¯tH production threshold occurs at s0 = (2mt +mH)2, alongside the particle-mass dependence of the parametrisation.Source: §A.1 · p. 42
• The scan variable β² vanishes at production threshold and approaches 1 in the high-energy limit, providing a convenient partonic-energy coordinate; see Ref..Source: §A.1 · p. 42