PaperlayerView exact revision on arXiv

Source-linked AI summary

A GPU compatible quasi-Monte Carlo integrator interfaced to pySecDec

S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk

arXiv:1811.11720v2hep-ph

TL;DR

Growing loop count and kinematic scales make multi-loop calculations increasingly complex, motivating numerical alternatives. The paper presents a CPU/GPU rank-1 shifted-lattice QMC integrator for stand-alone and pySecDec use; its authors report that the rules can outperform Monte Carlo and that GPUs often provide order-of-magnitude-or-more speedups.

  • Problem

    Multi-loop calculations become more complex as loop count and kinematic scales increase, while analytic results beyond multiple polylogarithms enter largely unexplored function classes.

    Source: §1 · p. 3

  • Method

    The paper provides a rank-1 shifted-lattice QMC integrator usable on CPUs and GPUs, standalone or within pySecDec, and combines QMC integration with importance sampling.

    Source: §1 · p. 4

  • Results

    Rank-1 shifted lattice rules can considerably outperform Monte Carlo, and GPUs often yield speedups of an order of magnitude or more.

    Source: §6 · p. 33

  • Takeaways & Limitations

    The implementation is demonstrated on integrals ranging from a two-loop pentagon to a six-loop two-point function.

    Source: §6 · p. 33

    Periodization can increase integrand variance exponentially with dimension, and larger lattice sizes do not always reduce integration error.

    Source: §2.3 · p. 10, §5.1 · p. 28

Abstract

from arXiv · show

The purely numerical evaluation of multi-loop integrals and amplitudes can be a viable alternative to analytic approaches, in particular in the presence of several mass scales, provided sufficient accuracy can be achieved in an acceptable amount of time. For many multi-loop integrals, the fraction of time required to perform the numerical integration is significant and it is therefore beneficial to have efficient and well-implemented numerical integration methods. With this goal in mind, we present a new stand-alone integrator based on the use of (quasi-Monte Carlo) rank-1 shifted lattice rules. For integrals with high variance we also implement a variance reduction algorithm based on fitting a smooth function to the inverse cumulative distribution function of the integrand dimension-by-dimension. Additionally, the new integrator is interfaced to pySecDec to allow the straightforward evaluation of multi-loop integrals and dimensionally regulated parameter integrals. In order to make use of recent advances in parallel computing hardware, our integrator can be used both on CPUs and CUDA compatible GPUs where available.

PROGRAM SUMMARY

pySecDec extracts singularities from dimensionally regulated integrals through sector decomposition, then numerically evaluates the resulting finite parameter integrals; the program also provides a QMC integrator and supports CPU/GPU execution.

Program summary

  • Sector decomposition factors singularities and yields finite unit-hypercube integrals for numerical evaluation of regulator-expansion coefficients.Integrable singularities can be handled by automatically choosing a suitable contour in the complex plane.Source: PROGRAM SUMMARY · p. 2
  • The QMC integrator is available as a stand-alone library and within pySecDec, with CUDA-compatible GPU support requiring Compute Capability 3.0 or greater.The program also supports CPU execution; running time ranges from seconds to days depending on problem complexity.Source: PROGRAM SUMMARY · p. 2
  • The listed external libraries include CUBA, FORM, and Normaliz, although a program mode does not require Normaliz.Other listed dependencies include catch, gsl, numpy, sympy, and Nauty.Source: PROGRAM SUMMARY · pp. 2–3

1 Introduction

The introduction motivates numerical integration as an alternative for increasingly complex multi-loop calculations, particularly when analytic results extend beyond multiple polylogarithms. It presents a CPU/GPU-capable QMC integrator, an importance-sampling method, and integration with pySecDec.

1 Introduction

  • Analytic calculations grow more challenging as loop count and kinematic scales increase, especially beyond multiple polylogarithms, where results may involve elliptic or hyper-elliptic functions.Numerical approaches need to be fast and accurate enough to provide results where analytic approaches reach their limits.Source: §1 · p. 3
  • Monte Carlo integration has error scaling of 1/√n, whereas QMC can scale as 1/n or better when its integrand requirements are met.GPU performance for sector-decomposition functions and applications to multi-scale integrals had already been demonstrated in prior work.Source: §1 · p. 4
  • The paper presents a stand-alone QMC integrator using multiple CPU cores and GPUs, adds it to pySecDec, and introduces importance sampling.The QMC implementation can also be used outside pySecDec.Source: §1 · p. 4

2 Description of the QMC method

The QMC method uses randomly shifted rank-1 lattice rules to estimate integrals and their error, with generating vectors, periodization, and variance reduction addressing practical convergence needs.

2.1 Quasi-Monte Carlo integration

  • For suitable weights, rank-1 shifted lattice rules can approach O(n^-1) convergence independently of dimension for Sobolev functions, while Korobov-space rates can approach O(n^-α).Random shifts permit an unbiased mean-square error estimate; typical applications require 10–20 shifts for a reliable estimate.Source: §2.1 · pp. 5–8

2.2 Generating vectors

  • Component-by-component construction selects each new generating-vector component to minimize worst-case error, at cost O(d n log n).The library provides generating vectors for fixed lattice sizes, constructed using equal product weights γ_i = 1/d.Source: §2.2 · p. 8

2.3 Transformations

  • Periodizing transforms improve smoothness for non-periodic integrands, but their variance increase can grow exponentially with dimension and spoil QMC convergence.For relatively low dimensions (d ≲ 8), higher-weight transforms may have variance costs counteracted by improved smoothness; the baker’s transform can be useful at moderate dimensions because it typically raises variance less than Korobov and Sidi transforms.Source: §2.3 · pp. 9–10

2.4 Variance reduction

  • The variance-reduction method fits a smooth transformation to the inverse CDF dimension-by-dimension, approximating the variance-reducing change of variables.The transformation parameters are fitted to a lattice-based CDF estimate; the authors find the ansatz works well for typical sector-decomposition functions, while other integrands may benefit from a modified fit.Source: §2.4 · pp. 10–11

3 Stand-alone usage of the integrator library

The stand-alone library can be installed as a C++11 header-only package and used to integrate functions with either accuracy-driven or speed-driven settings. The section also describes GPU requirements and a work-queue design for CPU/GPU execution.

3.1 Installation

  • The library is installed by including its C++11 header and linking against GSL; optional GPU support requires a CUDA compiler and compatible hardware with Compute Capability 3.0 or greater.The passage specifies GSL version 2.5 or greater.Source: §3.1 · pp. 11–12

3.2 Minimal example

  • A minimal program integrates the three-dimensional function f(x0, x1, x2) = x0x1x2 using default settings and prints the integral and its uncertainty.The example instantiates the integrator with a Korobov transform of weight 3; MAXVAR must be at least the largest number of integration parameters passed to that instance.Source: §3.2 · pp. 12–13

3.3 Usage

  • The library supports two use cases: meeting specified relative or absolute accuracy with little user input, or obtaining a result as quickly as possible when QMC performance is already understood.Source: §3.3 · p. 13

3.3.1 Case 1

  • For an accuracy-driven run, users set epsrel and epsabs and use maxeval to stop integration if the target error is not reached.The integrator chooses a periodizing transform, optionally fits a function for variance reduction, and repeats shifted-lattice estimates with larger lattices until the error goal or maxeval is reached.Source: §3.3.1 · pp. 14–15

3.3.2 Case 2

  • For a speed-driven run, users set minn and maxeval=1; the library evaluates minm shifted copies of the smallest available lattice with at least minn points and returns the result.If the result is unsatisfactory, the user can increase minn and retry.Source: §3.3.2 · p. 15
  • Similar-sized lattices can yield integral estimates with errors differing by an order of magnitude or more, making QMC scaling difficult to estimate straightforwardly.The authors suggest evaluating several lattice sizes when investigating scaling behavior.Source: §3.3.2 · pp. 15–16

3.3.3 Usage on GPUs

  • GPU use requires an integrand functor compatible with the target device, including supported C++ features, no host-only data access, and a host device call operator.GPU startup can add 1–3 seconds per device when the driver is not already initialized; virtual-architecture compilation may also become a dominant runtime fraction.Source: §3.3.3 · pp. 16–17

3.4 Design and Implementation

  • The central work queue balances work across CPU and GPU workers with different performance, but atomic locking can significantly hurt performance for quick integrands on many cores or devices.Balancing also depends on work packages being small enough that a slower worker does not finish much later than the others.Source: §3.4 · p. 17

4 Usage of the integrator library within pySecDec

This section explains how to install and use pySecDec’s QMC integrator on CPUs and CUDA-compatible GPUs, then illustrates it on multi-loop integrals with varied masses and kinematics. The examples include cases where fitted functions improve convergence and where QMC is practical for integrals that are difficult to compute analytically.

4.1 Installation

  • Install pySecDec after installing numpy and sympy; its make command also builds dependencies, including the QMC integrator library.The documentation points users to the online installation notes and getting-started guide.Source: §4.1 · p. 18

4.2.1 Using CPUs only

  • The CPU workflow defines the integral and options in Python, generates and builds the sector libraries, then runs numerical integration through the Python interface.Choose the QMC integrator in the integration script.Source: §4.2.1 · p. 19

4.2.2 Using GPUs and CPUs

  • CUDA compilation enables GPU use, and selecting the QMC integrator can use all CPU cores as well as all available GPUs.The GPU build uses Nvidia’s nvcc; specifying a virtual architecture can add just-in-time compilation time.Source: §4.2.2 · p. 19, §4.3.1 · p. 20

4.3.2 3-mass banana graph

  • The three-mass banana example evaluates a hyperelliptic integral, which the authors report can be evaluated as fast as integrals more accessible analytically.Its imaginary pole coefficients are numerically zero, with their accuracy limited by proximity to machine precision.Source: §4.3.2 · pp. 21–22

4.3.3 Non-planar 4-point function with massive propagators and massive legs of different mass / 4.3.4 Pentabox

  • The physical-region non-planar example reports a two-loop result with a massive loop and two differently massive external legs; the finite d = 6 − 2ϵ pentabox benefits considerably from a fit function.The pentabox’s finite dimensional version is presented as a more suitable numerical target than the 4 − 2ϵ version.Source: §4.3.3 · p. 22, §4.3.4 · pp. 22–23

4.3.5 Elliptic 2-loop integral / 4.3.6 Hyperelliptic 2-loop integral

  • For an elliptic physical-region integral, the authors state that QMC obtains correct digits that Monte Carlo cannot reach in a reasonable time; a massive hyperelliptic two-loop case is described as easily accessible numerically.The elliptic integral enters NLO Higgs-plus-jet corrections.Source: §4.3.5 · p. 23, §4.3.6 · p. 24

4.3.7 4-loop form factor example

  • For the 11-dimensional four-loop form factor, the Baker transform was needed for approximate 1/n scaling and outperformed the Korobov transform by not increasing integrand variance.The example also demonstrates sector decomposition for an integral with many loops and internal propagators.Source: §4.3.7 · pp. 24–25

4.3.8 2-loop Nbox / 4.3.9 6-loop bubble

  • The Nbox examples lack a Euclidean region, while symmetry finding reduces bubble6L sectors from more than 14,000 to 8,774 and requires geometric decomposition.Nbox split a and b use split=True, whereas split c is quasi-finite and does not need it; the iterative bubble decomposition leads to infinite recursion.Source: §4.3.8 · p. 25, §4.3.9 · p. 26

5 Profiling

The profiling tests show that QMC error scaling depends on the transform, lattice size, and integrand. Timing comparisons show QMC’s strengths and limits across test functions and loop integrals.

5.1 Scaling behaviour

  • For the O(ε4) contribution of a 3-loop massless form-factor integral, Korobov weight α = 1 gives relative precision 10^-9 at n ≈ 10^9 with approximate O(n^-1) scaling; α = 3 reaches 10^-14 at n ≈ 10^9 with approximate O(n^-2) scaling.The expected O(n^-3) asymptotic scaling was not observed for lattices with n ≈ 10^8.Source: §5.1 · p. 27
  • Larger n does not ensure smaller error: unusually uncertain lattice sizes occur for nearly all integrals, but which sizes are problematic depends on the integrand.Source: §5.1 · p. 28
  • For n ≲ 10^6, the Higgs+jet sector’s error scaled only as n^-1/2 and exceeded the true result; larger n gave O(n^-1) scaling and better than 0.1% precision.Importance sampling reduced error by about a factor of 3 at small n and more than 10 at large n; an integration-by-parts basis with n ≈ 10^6 sufficed for accurate results.Source: §5.1 · pp. 28–29

5.2 Timings for test functions

  • With εrel = 10^-8 and a ≲360-second limit, Divonne usually did not reach its second phase and often returned 2 correct digits; εrel = 10^-5 gave 4–5 digits in around 20 seconds.Source: §5.2 · p. 30
  • QMC often beat other Cuba integrators in digits and time, while GPU sampling took around 1 second versus 30–80 seconds on CPU; these simple tests poorly resemble sector-decomposed integrands.Source: §5.2 · pp. 30–31, Table 1 · p. 31
  • In d = 10, Baker-transform QMC could exceed the digits obtained in d = 8 with weight-3 Korobov; that Korobov transform can increase variance.Source: §5.2 · p. 32

5.3 Timings for loop integrals

  • QMC is generally faster than Vegas for obtaining many significant digits, especially when GPUs are available; Vegas can sometimes be faster for low-precision results.Source: §5.3 · p. 33

6 Conclusions

The paper presents a CPU- and GPU-capable quasi-Monte Carlo integrator for standalone use or through pySecDec, with importance sampling and examples spanning several multiloop integral types. Rank-1 shifted lattice rules often outperform Monte Carlo, and GPUs can provide order-of-magnitude-or-greater speed-ups in many cases.

6 Conclusions

  • The integrator uses rank-1 shifted lattice rules and runs on CPUs or GPUs, either as a standalone library or within pySecDec.The paper also describes a novel approach combining QMC integration with importance sampling.Source: §6 · p. 33
  • pySecDec examples cover a 2-loop pentagon, elliptic or hyperelliptic integrals, a 4-loop form factor, and a 6-loop 2-point function.The new pySecDec version also includes an improved algorithm for detecting sector symmetries.Source: §6 · p. 33
  • Rank-1 shifted lattice rules can considerably outperform Monte Carlo, and GPUs can provide speed-ups of an order of magnitude or more in many cases.The authors report that their implementation often computes more accurate digits in reasonable time than Vegas-like Monte Carlo; sector-decomposition integrands are typically smooth enough for O(1/n) scaling, unlike some NNLO phase-space integrands.Source: §6 · pp. 33–34
  • The authors suggest the publicly available implementation could automate, speed up, and improve accuracy for numerical multiloop amplitudes with several mass scales.The standalone integrator and pySecDec version are publicly available through the stated project repositories and documentation.Source: §6 · p. 34

A API documentation

The API documentation describes the QMC class's template parameters, variance-reduction fitting option, and support for compatible floating-point types.

A API documentation

  • The QMC class templates specify the integral's return and argument types, maximum variable count, and an optional integrand transform.The class also has template parameters for a fitting function, pseudo-random engine, and uniform real distribution.Source: A API documentation · pp. 34–35
  • The fitting-function parameter can fit the integrand's inverse cumulative distribution function dimension-by-dimension to reduce variance.Its default is fitfunctions::None::template type.Source: A API documentation · p. 35
  • Alternative floating-point types are supported when required standard functions are compatible or have suitable overloads.The documentation lists sqrt, abs, modf, pow, std::max, and std::min; non-real or non-complex types may also need error-calculation overloads.Source: A API documentation · p. 35
  • A Boost multiprecision binary-float example demonstrates instantiating the QMC class with a non-standard type.The example requires the Boost library.Source: A API documentation · p. 35

A.1 Public fields

The public fields configure logging, random-number generation, integration tolerances and limits, parallel execution, and variance-reduction fitting.

A.1 Public fields

  • The random generator can be seeded, and setting one CPU thread with devices={-1} disables multithreading for reproducibility.The library can direct output to a chosen file through its logger; verbosity ranges from 0 to 3.Source: §A.1 · pp. 35–36
  • The integrator's lattice size, number of random shifts, and maximum evaluations set integration workload and error estimation.Defaults include a minimum lattice size of 8191, 32 shifts, and 1,000,000 maximum evaluations; unavailable requested lattice sizes may be rounded up to the next available size.Source: §A.1 · p. 36
  • The error settings control absolute and relative goals; for complex results, errormode selects separate component goals or a goal based on the largest component.For real results, the target is |E-I| <= max(epsabs, epsrel*I); the default epsabs is 1e-7, and the default complex errormode is all.Source: §A.1 · pp. 36–37
  • Execution fields configure CPU threads, CUDA blocks and threads per block, and the CPU or GPU devices used.When GPUs are used, one additional CPU thread per device is launched for device communication; device id -1 represents CPUs.Source: §A.1 · p. 37
  • The generating-vector map determines available lattices; the default vectors support integrands with up to 100 dimensions.The documentation also describes an evaluateminn setting for points used to evaluate the integrand and, when enabled, fit the inverse cumulative distribution function.Source: §A.1 · pp. 37–38
  • Variance-reduction fit settings control how many points enter the fit and provide nonlinear least-squares tolerances and iteration parameters.fitstepsize selects every xth sorted point; its default is 10, while fitparametersgsl defaults to {}.Source: §A.1 · p. 38

gsl multifit nlinear parameters fitparametersgsl

The fitparametersgsl field uses parameters documented for GSL nonlinear least-squares fitting and defaults to an empty setting.

gsl multifit nlinear parameters fitparametersgsl

A.2 Public Member Functions

The public functions integrate an integrand, evaluate it on shifted lattice points, and fit a variance-reduced integrand. Their interfaces return integration results, sample data, or a fitted integrand as appropriate.

A.2 Public Member Functions

  • integrate evaluates func in d dimensions using the selected transform and returns the integral, estimated absolute error, lattice and shift counts, iterations, and evaluations.The result also records the largest lattice size and its number of shifts.Source: §A.2 · p. 39
  • evaluate samples func on a lattice of size at least evaluateminn and returns its generating vector, random shift, integrand values, lattice size, and sample-coordinate accessor.The accessor returns the argument used for a specified sample and integration variable.Source: §A.2 · p. 39
  • fit approximates the integrand's inverse cumulative distribution dimension-by-dimension and returns a functor representing the resulting integrand.The function being fitted is specified by the integrator's type F.Source: §A.2 · pp. 39–40
  • The integrand functor must expose its number of integration variables as a public member variable.The functions also call get next n where indicated by their interface descriptions.Source: §A.2 · pp. 39–40

A.3 Generating vectors

The program distributes rank-1 lattice generating vectors for specified lattice sizes and maximal dimensions, with custom vectors available when the supplied choices are insufficient.

A.3 Generating vectors

  • Distributed generating vectors cover maximal dimensions s ≤6 or s ≤100, and their types are summarized in Table A.1.The vectors use a Component-By-Component construction, computed partly with fastrank1pt.m and, for very large lattice sizes, the authors' FFTW-based CBC tool [64].Source: §A.3 · p. 40, Table A.1 · p. 40
  • The distributed vectors target Korobov spaces with smoothness α = 2, using the stated kernel, weights γ_i = 1/s, and parameters β_i = 1.The kernel is ω(x) = 2π2(x2 −x + 1/6), with i = 1, . . . , s.Source: §A.3 · p. 40
  • The generating vectors used by QMC are selected through the integrator's generatingvectors member variable.The example assigns a distributed CBC vector set to that member.Source: §A.3 · p. 40
  • Users can supply custom generating vectors when the available options, including up to 100 dimensions or 15173222401 lattice points, are insufficient.The passage illustrates setting a map with vectors z = (1, 3) for n = 7 and z = (1, 7) for n = 11.Source: §A.3 · pp. 40–41

A.4 Integral Transforms

The integrator lets users select both periodizing transformations and fit functions. The transformations apply weights before integration, while the fit functions approximate the integrand's inverse cumulative distribution dimension-by-dimension to reduce variance.

A.4 Integral Transforms

  • A Korobov transform with r0 = 5 and r1 = 3 applies a weight to the integrand before integration.Available periodizing transformations are listed in Table A.2.Source: §A.4 · p. 41, Table A.2 · p. 41

A.5 Fit Functions

  • Fit functions approximate the integrand's inverse cumulative distribution dimension-by-dimension, and a PolySingular fit function reduces variance before integration.Possible fit-function types are listed in Table A.3.Source: §A.5 · pp. 41–42, Table A.3 · p. 41
Loading 1811.11720v2…