Integration-by-parts reductions of Feynman integrals using Singular and GPI-Space

TL;DR

This paper presents a novel algebro-geometric IBP reduction method for complex Feynman integrals, leveraging parallel computation frameworks to automate and optimize the reduction process.

Contribution

It introduces a new IBP reduction approach combining algebraic geometry, sparse linear algebra, and parallel computing with Singular and GPI-Space.

Findings

Successfully reduced two-loop five-point nonplanar double-pentagon integrals.

Automated and parallelized the IBP reduction process.

Potential for simplifying IBP coefficients in transcendental bases.

Abstract

We introduce an algebro-geometrically motived integration-by-parts (IBP) reduction method for multi-loop and multi-scale Feynman integrals, using a framework for massively parallel computations in computer algebra. This framework combines the computer algebra system Singular with the workflow management system GPI-Space, which is being developed at the Fraunhofer Institute for Industrial Mathematics (ITWM). In our approach, the IBP relations are first trimmed by modern algebraic geometry tools and then solved by sparse linear algebra and our new interpolation methods. These steps are efficiently automatized and automatically parallelized by modeling the algorithm in GPI-Space using the language of Petri-nets. We demonstrate the potential of our method at the nontrivial example of reducing two-loop five-point nonplanar double-pentagon integrals. We also use GPI-Space to convert the basis of IBP reductions, and discuss the possible simplification of IBP coefficients in a uniformly transcendental basis.

Tables

Table 1. Table 1: Timings and efficiency for the cut { 1 , 3 , 4 , 5 } 1 3 4 5 \{1,3,4,5\} . We use the same algorithm for all core counts. The single core run serves as a reference.

			relative
nodes	cores	runtime	speedup	efficiency
1	1	122857.6	1.201	1.201
1	15	9837.8	15.000	1.000
2	30	4954.8	29.7822	0.992
4	60	2625.4	56.2058	0.936
8	120	1341.3	110.014	0.916
14	210	952.3	154.958	0.737
15	225	705.6	209.132	0.929
16	240	694.3	212.514	0.885
29	435	611.8	241.199	0.554
30	450	385.4	382.856	0.850
32	480	379.9	388.336	0.809
40	600	367.7	401.307	0.668
48	720	363.2	406.195	0.564

Table 2. Table 2: The IBP relations generated on each cut by the module intersection method. We used finite-field methods to pick linearly independent and necessary IBP relations to reduce all target integrals. The size is the output file size on disk before reduction. The numbers d 2 subscript 𝑑 2 d_{2} and d 3 subscript 𝑑 3 d_{3} are the maximal degrees in the reduced IBP relations for c 2 subscript 𝑐 2 c_{2} and c 3 subscript 𝑐 3 c_{3} , respectively.

Cut	# relations	# integrals	size	$d_2$	$d_3$
{1,5,7}	1134	1182	0.77 MB	21	22
{1,5,8}	1141	1192	0.85 MB	18	18
{1,6,8}	1203	1205	1.1 MB	19	30
{2,4,8}	1245	1247	1.1 MB	35	24
{2,5,7}	1164	1211	0.84 MB	26	18
{2,6,7}	1147	1206	0.62 MB	16	17
{2,6,8}	1126	1177	0.83 MB	16	18
{3,4,7}	1172	1221	0.78 MB	19	18
{3,4,8}	1180	1226	1.0 MB	19	22
{3,6,7}	1115	1165	0.82 MB	21	28
{1,3,4,5}	721	762	0.43MB	14	14

Equations

$$0=\int \frac{{\mathrm{d}}^D{\ell }_1}{\mathrm{i}{\pi }^{D/2}}\dots \frac{{\mathrm{d}}^D{\ell }_L}{\mathrm{i}{\pi }^{D/2}}\sum _{j=1}^L\frac{\partial }{\partial {\ell }_j^{\mu }}\frac{v_j^{\mu }}{D_1^{{\nu }_1}\cdots D_m^{{\nu }_m}},$$

$$I[n_1,\dots ,n_m]=\int \frac{{\mathrm{d}}^D{l}_1}{i{\pi }^{D/2}}\dots \frac{{\mathrm{d}}^D{l}_L}{i{\pi }^{D/2}}\frac{1}{D_1^{n_1}\cdots D_m^{n_m}},$$

$$I[n_1,\dots ,n_m]$$

$$P=\det G\left(\begin{array}[]{cccccc}l_{1},&\ldots&l_{L},&p_{1},&\ldots&p_{E}\\ l_{1},&\ldots&l_{L},&p_{1},&\ldots&p_{E}\end{array}\right)\,.$$

$$U=\det G\left(\begin{array}[]{ccc}p_{1},&\ldots&p_{E}\\ p_{1},&\ldots&p_{E}\end{array}\right),\quad C_{E}^{L}=J\frac{\pi^{\frac{L-m}{2}}}{\Gamma(\frac{D-E-L+1}{2})\ldots\Gamma(\frac{D-E}{2})}\,,$$

$$0$$

$$\displaystyle\bigg{(}\sum_{i=1}^{m}a_{i}(z)\frac{\partial P}{\partial z_{i}}\bigg{)}+b(z)P=0\,,$$

$$a_i(z)=b_i(z)z_i,i=1,\dots ,\kappa ,$$

$$M_1\cap M_2.$$

$$c_i\mapsto a_i,1\le i\le k_1,$$

$$c_j\mapsto b_j,b_j\in \mathbb{Z},k_1<j\le k,$$

$$(d_1+1)\times (d_2+1)\times \dots \times (d_{k_1}+1),$$

$$(M{|}_{c_1\mapsto a_1,\dots ,c_r\mapsto a_r}{)}_{\text{red}}.$$

$$M_{\text{red}}{|}_{c_1\mapsto a_1,\dots ,c_r\mapsto a_r}$$

$$M_{\mathrm{red}}{|}_{c_1\mapsto a_1^{(0)}},\dots ,M_{\mathrm{red}}{|}_{c_1\mapsto a_1^{(d_1)}},$$

$$M_{\text{ref}}:={(M{|}_{c_{r+1}\mapsto q_{r+1},\dots ,c_k\mapsto q_k})}_{\text{red}}.$$

$$M^{(j)}:=M{|}_{c_1\mapsto p_1^{(j)},\dots ,c_{j-1}\mapsto p_{j-1}^{(j)},c_{j+1}\mapsto p_{j+1}^{(j)},\dots ,c_k\mapsto p_k^{(j)}}$$

$$(M_{c_1\mapsto p_1^{\prime },\dots ,c_m\mapsto p_r^{\prime }}){|}_{\text{red}}.$$

$$\lceil \frac{450}{\text{number of CPUs}}\rceil .$$

$$T_{\text{ideal}}(n)=(1-f)\cdot T(1)+\frac{f\cdot T(1)}{n},$$

$$S_{\text{ideal}}(n)=\frac{n}{(1-f)\cdot n+f},E_{\text{ideal}}(n)=\frac{1}{(1-f)\cdot n+f}.$$

$$dI=d{x}_1\wedge d{x}_2\wedge ...\wedge d{x}_{n}R(x_1,...,x_n),$$

$$dI=\sum _i{c}_{i}d\text{log}{f}_{i,1}\wedge ...\wedge d\text{log}{f}_{i,n},$$

$$\frac{dx}{x-a}\wedge \Omega =d\log (x-a)\wedge \Omega ,$$

$$F=AI,$$

$$\tilde{I}=TI.$$

$$D_1$$

$$D_5$$

$$D_9$$

$$I[1,1,1,1,1,1,1,1,0,0,-4],I[1,1,1,1,1,1,1,1,0,-1,-3],I[1,1,1,1,1,1,1,1,0,-2,-2],$$

$$I[1,1,1,1,1,1,1,1,0,-3,-1],I[1,1,1,1,1,1,1,1,0,-4,0],I[1,1,1,1,1,1,1,1,-1,0,-3],$$

$$I[1,1,1,1,1,1,1,1,-1,-1,-2],I[1,1,1,1,1,1,1,1,-1,-2,-1],I[1,1,1,1,1,1,1,1,-1,-3,0],$$

$$I[1,1,1,1,1,1,1,1,-2,0,-2],I[1,1,1,1,1,1,1,1,-2,-1,-1],I[1,1,1,1,1,1,1,1,-2,-2,0],$$

Full text

[TABLE]

1 Introduction

With the success of Large Hadron Collider (LHC) Run II and the upcoming LHC run III, high precision background computation, especially next-to-next-to-leading-order (NNLO) scattering computation, is crucial for the interpretation of experimental results. In recent years, great progress has been made in multi-loop scattering amplitude calculations, for instance, in the case of $2\rightarrow 3$ processes Badger:2013gxa ; Gehrmann:2015bfy ; Badger:2017jhb ; Abreu:2017hqn ; Abreu:2018aqd ; Abreu:2018jgq ; Boels:2018nrr ; Gehrmann:2018yef ; Badger:2018enw ; Abreu:2018zmy ; Chicherin:2018yne ; Chicherin:2019xeg ; Abreu:2019rpt ; Abreu:2019odu ; Hartanto:2019uvl . The progress is due to modern developments of scattering amplitudes, like the integrand construction method Zhang:2012ce ; Mastrolia:2012an , canonical integrals Henn:2013pwa ; Henn:2014qga , numeric unitarity Ita:2015tya ; Abreu:2017xsl , bootstrap methods Dixon:2011pw ; Dixon:2013eka ; Dixon:2014iba ; Caron-Huot:2016owq ; Dixon:2015iva ; Dixon:2016nkn ; Chicherin:2017dob ; Caron-Huot:2019vjl , reconstruction using finite fields vonManteuffel:2014ixa ; Peraro:2016wsq ; Klappert:2019emp ; Peraro:2019svx and new ideas in the integration-by-parts (IBP) reduction. The latter is the main topic to be discussed in this paper.

Frequently, when computing scattering amplitudes, IBP reduction is a crucial and bottleneck step. It is a fundamental tool for both the reduction of integrals to master integrals (MIs), and for computing the master integrals themselves using the differential equation method. IBP relations (IBPs) are derived from integrating a total derivative CHETYRKIN1981159 ,

[TABLE]

where the $v_{i}^{\mu}$ are polynomials in the loop momenta $\ell_{i}$, the $D_{i}$ are the inverse propagators, and $D$ is the spacetime dimension.

The standard approach to obtain IBP reductions, by which we are able to express an integral as a linear combination of a finite number of MIs, is to generate sufficiently many IBP relations, and then use the Laporta algorithm Laporta:2001dd to solve the associated linear system. The algorithm works by imposing an ordering on the different integral families and solving recursively. There exist multiple public and private implementations of this approach Smirnov:2008iw ; Smirnov:2014hma ; Smirnov:2019qkx ; Maierhoefer:2017hyi ; Maierhofer:2018gpa ; vonManteuffel:2012np ; Klappert:2019emp , which usually generates a large linear system to be solved.

In the case of a system of IBPs which does not have double propagators Gluza:2010ws ; Schabinger:2011dz ; Larsen:2015ped , however, we obtain a much smaller linear system. The IBPs without double propagators are physically related to dual conformal symmetries Bern:2017gdk . A significant simplification can be made by using unitarity methods, where by considering a spanning set of generating cuts it is possible to reduce the size of the IBP system. This requires prior knowledge of a basis of MIs. Such a basis can be obtained by running the Laporta algorithm with constant kinematics, or by using specialized programs such as Mint Lee:2013hzt or Azurite Georgoudis:2016wff . (Note that the dimension of a basis of integrals can also be obtained by studying the parametric annihilators Bitoun:2017nre .) There is also the important technique Chawdhry:2018awn of simultaneously nullifying all master integrals except one, which often makes large-scale linear reductions feasable.

Besides the advances in purely analytical methods in recent years, there has been a lot of work towards numerical implementations of the generation of IBPs. The idea is to utilize either integer values or finite-field values for the kinematical invariants vonManteuffel:2014ixa ; Peraro:2016wsq ; Smirnov:2019qkx , depending on the difficulty of the problem, and then to run the same reduction several times for reconstruction. This method has been very successful in tackling difficult problems. Furthermore, it is possible to numerically generate and reduce the IBP relations, and, while skipping the IBP coefficient reconstruction, directly carry out an amplitude reconstruction. (For examples, see Badger:2018enw ; Abreu:2018zmy ; Abreu:2019rpt ; Badger:2019djh ). In this paper, we in particular present our own implementation of a semi-numeric rational interpolation method, see Appendix A for more details.

Furthermore, new approaches were developed recently to obtain the reduction directly, without generating IBP relations from total derivatives. In Kosower:2018obg , the direct solution method was presented to derive recurrence relations for Feynman integrals and solve them analytically with arbitrary numerator degree. One very promising progress is based on the intersection theory of differential forms in the Baikov representation Mastrolia:2018uzb ; Frellesvig:2019kgj ; Frellesvig:2019uqt . This approach calculates the master integral coefficients from intersection numbers. There is also a very intuitive approach to reduce Feynman integrals by considering the $\eta$ expansion of the Feynman prescription Liu:2018dmc ; Liu:2017jxz ; Wang:2019mnn . Using this approach, the scaling of the reduction computation depends only linearly on the number of master integrals. Furthermore, it is possible to determine two-loop planar diagram IBP coefficients directly from the Baikov representation Kardos:2018uzy .

In this paper, we present our new powerful IBP reduction method based on:

Computational algebraic geometry. We apply the module intersection method from Boehm:2018fpv ; Boehm:2017wjc , modified by using a suitably chosen degree bound for the Gröbner basis computation, to efficiently generate a small IBP system, without double propagators (or IBPs with a given bound on the propagator exponents). 2. 2.

A modern framework of massively parallel computations in computer algebra which combines the computer algebra system Singular singular with the workflow management system GPI-Space GPI . We have completely automatized our approach and make our algorithms run automatically in parallel on high performance computing clusters. In this way, IBP results can be obtained in an efficient, reliable and scalable way. Our implementation can automatically determine the minimal number of points needed for interpolating the IBP coefficients, it can identify possible “bad” points, add more points, if necessary, and interpolate the final result.

We demonstrate the power of our method by reducing the two-loop five-point nonplanar double pentagon diagram analytically, up to numerator degree $4$. This is a nontrivial test since the diagram has a complicated topology and there are five symbolic Mandelstam variables as well as the spacetime variable $D$.

Furthermore, we start to look at the possible simplification of IBP coefficients by converting the master integral basis. In this paper, we test the conversion to a “dlog” basis WasserMSc , a special case of the canonical basis Henn:2014qga . We find that for the double pentagon diagram above, the size of the IBP coefficients reduces significantly from the byte size $\sim 2.0$G in the Laporta basis to $\sim 0.48$G in the dlog basis on disk, that is, by $76\%$. The master integral basis conversion computation is also automated by the Singular-GPI-Space framework.

Our paper is structured as follows. In Section 2, we present the general background on how to generate simple and trimmed IBP systems using computational algebraic geometry and finite-field methods, as well as the improvement on the algorithm in Boehm:2018fpv . In Section 3, we give a short overview on how we use Singular in conjunction with GPI-Space. In Section 4, we describe how to model our algorithm in the Singular-GPI-Space framework, and discuss timings and scaling of the algorithm, focusing on the double pentagon diagram. This, in particular, demonstrates the potential of the Singular-GPI-Space framework for applications in high-energy physics. In Section 5, we review the algorithmic computation of a dlog basis which has uniform transcendental weight, and we comment on how to convert coefficients from the Laporta basis to the dlog basis. In Section 6, we study the working example of our implementation, the double pentagon graph, in detail. We discuss the analytic IBP reduction and the conversion of IBP coefficients to the dlog basis. Finally we present a summary and conclusion of this paper.

The result of our IBP reductions can be downloaded from the following links: Whereas

https://www.dropbox.com/s/1ubdhcyhe8e4pwy/IBPmatrix_Laporta_basis.tar.gz

provides the IBP coefficients in the Laporta basis with the scale $s_{12}=1$,

https://www.dropbox.com/s/e6t4evftkfo95pr/IBPmatrix_dlog_basis.tar.gz

contains the IBP coefficients in the dlog basis with the scale $s_{12}=1$.

For the convenience of the reader, we also present the IBP coefficients in the dlog basis with the full scale dependence:

https://www.dropbox.com/s/dnkr6h5t3vik2r0/IBPmatrix_dlog_basis_scaled.tar.gz

We encourage researchers in the high energy community to send us IBP reduction problems (mailto: alessandro.georgoudis@physics.uu.se) for cutting-edge precision calculations and the further sharpening of our new reduction method.

2 The module intersection method reloaded

In this section, we present a refined version of the approach of using module intersections to trim IBP systems. For the detailed account of the module intersection IBP reduction method, we refer to Boehm:2018fpv .

2.1 Module intersection

The Feynman integrals under consideration are labeled as

[TABLE]

where $L$ is the loop order and the $l_{i}$’s are the loop momenta. We have $E$ independent external vectors that we label as $p_{1},...,p_{E}$. We assume that the Feynman integrals have been reduced on the integrand level, and set $m=LE+L(L+1)/2$ which equals the number of scalar products in the configuration.

For us it is convenient to use the Baikov representation Baikov:1996rk ; Lee:2013hzt for IBP reductions,

[TABLE]

Here, $P$ is the Baikov polynomial, which can be written as a Gram determinant,

[TABLE]

Moreover, $U$ and $C_{E}^{L}$ are the Gram determinant respectively constant factor below:

[TABLE]

where $J$ is a constant Jacobian. The factors $U$ and $C_{E}^{L}$ are irrelevant for the IBP relations.

As in Lee:2014tja ; Ita:2015tya ; Larsen:2015ped , the IBP relations in the Baikov representation are of type

[TABLE]

where each $a_{i}(z)$ is a polynomial in the variables $z_{1},\dots,z_{m}$. Note that $P$ vanishes on the boundary of the Baikov integration domain, so this form of IBP identities does not have surface terms.

Suppose we wish to reduce an integral family with $n_{j}\leq 0$, $j=\kappa+1,\ldots,m$, for some $\kappa$. That is, we face integrals with the inverse propagator product $1/(D_{1}\ldots D_{\kappa})$ and the sub-topology integrals. We use the idea of restricting to IBP systems without double propagators Gluza:2010ws , choosing suitable $a_{i}(z)$ to prevent the appearance of double propagators in (2.5). In the Baikov representation, we also need to avoid total derivatives with dimension shifts Ita:2015tya ; Larsen:2015ped . These constraints translate into syzygy equations of the following type:

[TABLE]

where $b(z)$ and the $b_{i}(z)$ are also polynomials in $z_{i}$’s. Relation (2.6) avoids dimension shifts of the integrals, while (2.7) ensures that there is no double propagator for $D_{i}$ if the initial index $n_{i}=1$ in (2.5). The goal is to find such polynomials $a_{i}(z)$, $b(z)$, and $b_{i}(z)$. Since we require polynomial solutions, this is not a linear algebra problem, but a computational algebraic geometry problem.

We use the module intersection method from Zhang:2016kfo ; Boehm:2018fpv to solve (2.6) and (2.7) simultaneously. Note that the analytic generators of all solutions of (2.6) can be directly written down via either the canonical IBP vector method Ita:2015tya or the Gram matrix Laplace expansion method Boehm:2017wjc 111We learned the Laplace expansion method from Roman Lee, and proved its completeness via the Gulliksen-Negard/Jozefiak exact sequence in Boehm:2017wjc .. The relations in (2.7) can be trivially expressed as a module membership condition. Hence without any algorithmic computation, we know the individual solutions for (2.6) and (2.7), respectively. These form polynomial submodules $M_{1}$ respectively $M_{2}$ of $R^{m}$ over the polynomial ring $R=\mathbb{Q}(c_{1},\ldots,c_{k})[z_{1},\ldots,z_{m}]$ (where the variables $c_{1},\ldots c_{k}$ collect the Mandelstam variables and the mass parameters). The task is then to compute

[TABLE]

This module intersection can be obtained by computing a module Gröbner basis in a particular ordering Boehm:2018fpv . One decisive strategy is the localization technique, which allows us to compute $M_{1}\cap M_{2}$ over the polynomial ring $\tilde{R}=\mathbb{Q}[c_{1},\ldots,c_{k},z_{1},\ldots,z_{m}]$. In this manner, we treat kinematic variables in the same way as the Baikov variables. This greatly speeds up the intersection computation for multi-scale problems, but results in a redundant generating system. The latter can be trimmed further by importing the result back to $R^{m}$ and removing redundant generators by checking the leading monomials. This is powered by Singular’s command simplify. Once $M_{1}\cap M_{2}$ is obtained, we know all simultaneous solutions for (2.6) and (2.7), and can use (2.5) to get IBPs without double propagators.

We emphasize that, although (2.6) and (2.7) were originally designed for IBPs without double propagators, the solutions of (2.6) and (2.7) can be used to simplify IBP systems with double or multiple propagators. Using these solutions $a_{i}(z)$, the resulting IBP system does not introduce integrals with higher powers of propagators, and hence also greatly decreases the size of the IBP system.

Frequently, instead of computing IBPs directly, we compute IBPs on spanning cuts and assemble the full IBPs afterwards. This amounts to setting some of the $z_{i}$ to zero in (2.6) and (2.7). For details on IBPs on cuts using the Baikov representation, we refer to Boehm:2018fpv .

Compared to the approach in Boehm:2018fpv , we present the following new features of the module intersection method in this paper:

• •

When we compute the intersection $M_{1}\cap M_{2}$, instead of finding a full generating system, we heuristically impose a polynomial degree bound in the computation. Then we reduce the resulting IBPs over finite fields to test if we already have all the IBP relations needed. If the IBP relations are insufficient, we increase the degree bound and repeat the computation. This approach speeds up the intersection computation dramatically in many cases. In practice, we use the option degbound in the computer algebra software Singular to set the degree bound.

• •

In the approach of Boehm:2018fpv , the module intersection was only computed for the top sector, which, for the hexagon box diagram, turned out to be sufficient for reducing integrals to a master integral basis. However, in this paper, we compute the module intersection for the top sector and also all subsectors. This approach may, in general, generate more IBP relations. Via linear algebra trimming as discussed in the next subsection, this approach eventually gives a block triangular linear system and makes the linear reduction easier.

2.2 Linear reduction

For the simplified IBP system arising from the module-intersection method, we use our own linear reduction algorithm to reduce the IBP system. The steps are:

Trim the linear system in two stages: (a) Set all the kinematic variables to integer values, and use linear algebra over a finite field to find the independent IBP relations. (b) Again over a finite field, carry out the reduction. From the intermediate steps, determine a sufficient subset of IBP relations for reducing the target integrals. These operations are powered by the finite field computation tool SpaSM spasm . 2. 2.

Remove the overlap between two different cuts and simplify the linear system: If two cuts have a common master integral, use the idea from Chawdhry:2018awn to set the master integral to zero in the IBP system of one of the two cuts. This will later on dramatically simplify the IBP reduction for the cut. 3. 3.

For the linear system simplified by the first two steps, we use our own Singular row reduce echelon form (RREF) code over function fields to reduce the target integrals to master integrals. Our code applies both row and column swaps for finding the optimal pivots. Note that column swaps change the set of master integrals. After the RREF computation, we convert the new master integrals to the original master integrals. We have observed that this approach is in general much faster than fixing the column ordering and directly reducing the target integrals to the original master integrals.

For difficult IBP reduction computations, we use a “semi-numeric” approach: This approach sets several but usually not all of the kinematic variables for the reduction computation to numeric values (that is, to constant integers). Without loss of generality, for the kinematic variables $(c_{1},\ldots,c_{k})$, we set

[TABLE]

for some $k_{1}<k$ and some $a_{i}\in\mathbb{Z}$.

The actual degree of the coefficients in these variables can be decided by a univariate analytic computation (that is, we set all but one of the $c_{i}$ to constant values). For example, we may pick the dimension $D$ and all parameters $c_{i}$ except $c_{1}$ as random integers, and then carry out the reduction. This computation is much easier than the actual IBP reduction with fully analytic parameters. From the reduction, we determine the degree of $c_{1}$ in the final IBP reduction coefficients. Proceeding similarly for each $i$, we find the degree of each $c_{i}$. This determines the minimal number of semi-numeric points for the subsequent interpolation step. (See Peraro:2016wsq for an alternative way of finding the degree of each parameter in a rational function.)

After accumulating enough points, we collect the semi-numeric reduction results and interpolate to get the final IBP reduction coefficients. To do this, we first run step 3 above for a semi-numeric set of parameters, find the optimal pivots and record the row/column swap history as a trace of our computation. For other numeric values, we always use the same trace to ensure the relatively uniform running time of the computation.

In practice, we use our rational function interpolation algorithm described in Appendix A. We do a reduction computation, with a carefully chosen semi-numeric reference point,

[TABLE]

and $c_{1},\ldots c_{k_{1}}$ symbolic. Using the reference point result, we convert the rational function interpolation problem to individual polynomial interpolation problems for the numerators and denominators. With this approach, the number of “semi-numeric” computations is

[TABLE]

where the $d_{i}$, for $1\leq i\leq k_{1}$, are the maximal degrees of the $c_{i}$ in the numerator and denominator polynomials in the RREF matrix. This algorithm is also implemented in Singular.

For the semi-numerical reduction and interpolation, we need to parallelize our computations in an efficient way. Furthermore, with semi-numeric points, we may have some bad points in the reduction or interpolation. In order to make use of massively parallel computations in an efficient way, and to automize the workflow for the replacement of bad points, we use the modern workflow management system GPI-Space, in conjunction with the computer algebra system Singular. We will discuss the ideas behind this approach in the subsequent section.

3 Massively parallel computations using Singular and GPI space

Large scale calculations such as row reductions of IBP identities in the case of Feynman diagrams which are relevant to current research in high-energy physics, are only feasible by using parallel computing on high-performance clusters. The computer algebra methods applied in this context require to model algorithms which rely on sub-computations with time and memory requirements that are difficult to predict. This is due, for example, to the behaviour of Buchberger’s algorithm for finding Gröbner bases: Although this algorithm performs well in many practical examples of interest, its worst case complexity is doubly exponential in the number of variables Mayr1982TheCO . Nevertheless it turned out recently BDFPRR ; mptopcom that massively parallel methods, which have been a standard tool in numerical simulation for many years, can also be applied successfully in symbolic computation. Proposing the general use of massively parallel methods in computer algebra, we describe our ongoing effort in this direction which is based on connecting the computer algebra system Singular for polynomial calculations with the workflow management GPI-Space. The latter consists of a scheduler distributing the actual computations to workers in the cluster, a virtual memory layer to facilitate communication between the workers, and a workflow management system which relies on modeling algorithms in terms of Petri nets.

In its basic form, a Petri net is a directed bipartite graph with two kinds of nodes: While a place can hold a number of indistinguishable (structure-less) tokens, a transition may fire if each input place contains at least one token (we then say that the transition is enabled). When fired, a transition consumes one token from each input place and puts one token on each output place. See Figure 1 for an enabled transition and its firing, and Figure 2 for a transition which is not enabled. In the figures, places are shown as circles, transitions as rectangles, and tokens as black dots.

The execution of a Petri net is non-deterministic: At each step, a single random enabled transition is chosen to fire. We have observed that the randomized reformulation of deterministic algorithms in computer algebra in terms of Petri nets can lead to a more consistent and predictable behavior throughout the course of the computation.

In our approach, we model the coarse-grained structure of an algorithm in terms of a Petri net. The transitions call procedures from the C-library version of Singular to do the actual computations. The result of this setup is a flexible framework for massively parallel computations in computational algebraic geometry (similar setups are possible using C-libraries of computer algebra systems aiming at possibly different application areas). Our framework has, for example, already been used to implement a non-singularity test for algebraic varieties Ristau ; BDFPRR , the computation of combinatorial objects in geometric invariant theory Reinbold , and the computation of tropical varieties associated to algebraic varieties Bendle .

For the efficient use in practical programming, the basic concept of a Petri net has to be extended. Here, GPI-Space provides multiple additional features:

• •

Modeling complex algorithms just by the use of structure-less tokens is not very efficient. In GPI-Space, tokens can have a data type and hold actual data. In fact, it is often more efficient if the tokens just hold a reference to a storage place for the data (in memory or in the file system). Using the shared memory subsystem of GPI-Space or the powerful file systems of modern high-performance clusters, computations can then scale far beyond the limitations of a single machine.

• •

The firing of a transition may be subject to conditions which have to be fulfilled by the input tokens.

• •

Transitions in practice involve computations which take time. The properties of Petri nets allow us to execute different enabled transitions at the same time (task parallelism) and to execute multiple instances of the same transition in parallel, provided the input places hold multiple tokens (data parallelism). In Figure 3, the transitions $f_{1}$ and $f_{2}$ can fire in parallel, and, if the input place of $f_{i}$ holds multiple tokens, then $f_{i}$ can fire in multiple instances.

We have observed that some algorithms in computer algebra scale in a superlinear way when implemented in parallel as a Petri net. The reason is that then, at run time, the algorithms can automatically determine from a given set of paths a path which leads to the solution in the fastest possible way (see (BDFPRR, , Section 6.2)).

In the next section, we illustrate the use of the Singular-GPI-Space framework for applications in high-energy physics by modeling our IBP reduction algorithm.

4 Parallel matrix reduction as a Petri net

In this section, we first describe how to model the parallel IBP reduction algorithm in terms of a Petri net. Focusing on the cut $\{1,3,4,5\}$ of the two-loop five-point nonplanar double pentagon diagram, we then discuss timings and scaling of the algorithm to indicate the practical use and significant potential of the Singular-GPI-Space framework for algorithmic problems in high-energy physics.

4.1 General structure of the algorithm

Our approach includes a massively parallel execution of row-reductions over function fields, where a number of parameters has been replaced by integers, followed by a parallel interpolation step to reconstruct the dependency on these parameters.

So the task is to find the reduced row-echelon form $M_{\text{red}}$ of a large linear system of equations, given as a matrix $M$ over the rational function field $\mathbb{Q}(c_{1},\dots,c_{k})$. Since applying Gaussian elimination directly is not feasible, we instead proceed by substituting, say, the first $r$ parameters by the coordinates of a point $a\in\mathbb{Z}^{r}$, and then by computing the reduction

[TABLE]

We refer to Section 2.2 above for details on how we handle this reduction step. To determine the number of interpolation points required to reconstruct the dependency on $c_{1},\dots,c_{r}$, we find bounds for the degrees of numerators and denominators for each parameter by doing a univariate row reduction (that is, all but one of the parameters are set to be numeric). After the reduction, we check that the resulting matrix is equal to the desired result

[TABLE]

by normalizing it relative to a previously computed reference matrix with $c_{r+1},\dots,c_{k}$ constant, and performing degree checks using the exact degrees obtained from the univariate calculations. These steps are described in more detail in Appendix A. The final result $M_{\text{red}}$ is then found by iteratively combining the reduced matrices via univariate interpolation (see again Appendix A).

Let $d_{1},\dots,d_{r}$ be degree bounds for the entries of $M_{\text{red}}$ in the parameters $c_{1},\dots,c_{r}$, respectively. To obtain $M_{\text{red}}$ by interpolation, we need $d_{1}+1$ matrices over $\mathbb{Q}(c_{2},\dots,c_{k})$ of the form

[TABLE]

for $d_{1}+1$ values $a_{1}^{(0)},\dots,a_{1}^{(d_{1})}\in\mathbb{Z}$. Similarly, to obtain any one of the above matrices, we need $d_{2}+1$ matrices over $\mathbb{Q}(c_{3},\dots,c_{k})$. Continuing inductively, this process ends with matrices defined over $\mathbb{Q}(c_{r+1},\dots,c_{k})$, which are then computed by reduction with $c_{1},\dots,c_{r}$ numeric. This tree-like dependency structure is depicted in Figure 4.

4.2 Managing the interpolation

We model the current status of the interpolation process in a tree-like data structure corresponding to that from Figure 4, with references to the reduction results at the leaves, and references to the interpolation results at the other nodes. Within GPI-Space, reductions and interpolations are executed according to this data structure. The tree is generated as soon as the degree bounds $d_{1},\dots,d_{r}$ are known, and it is extended if the algorithm requires additional data points due to the occurrence of bad interpolation points.

4.3 Description of the Petri net

Figure 5 depicts the Petri net that implements the complete reduction algorithm. Going beyond the standard syntax introduced in Section 3, dashed arrows stand for read-only access, that is, the data in the respective places is not consumed. The dotted arrows illustrate read and write access to the interpolation tree described in Section 4.2. A transition can be annotated by conditions which indicate that the transition can only fire by consuming tokens for which the conditions evaluate to true.222When formulating conditions in the Petri net, we use the name of a place and a token on the place interchangeably. In the following, we describe the individual structures of the net:

Input token:

The net is initialized with one token: A token on the place $I$, which holds references to the following input data:

• •

The input linear relations, which are given as a matrix $M$ over the rational function field $\mathbb{Q}(c_{1},\dots,c_{k})$.

• •

The vector of indices of the parameters which will be interpolated (in the following we assume that these indices are $1,\dots,r$).

• •

The vector of indices of the target variables.

• •

Optionally: A precomputed trace for the reduction step (consistent with the targets). In the Petri net, the trace is referred to as $I$.trace (we use the usual dot-notation for sub-data structures). Note that the trace fixes the variables corresponding to the master integrals.

Transition trace:

If the token on $I$ does not contain a trace, then trace is enabled, computes a trace for the linear reduction (see Section 2.2) and returns a copy of $I$ with the trace included.

Transition copy:

If the token on $I$ already contains a trace, then copy is enabled and simply passes the token on $I$ through.

Transition init:

This transition takes the input token, which was produced by either trace or copy, and pushes it onto $I_{t}$. This way, the input data on $I_{t}$ is guaranteed to contain trace data. It additionally enables the transitions degrees and reference.

Transition reference:

This generates a random substitution point $q=(q_{r+1},\dots,q_{k})$ with values for all parameters which will not be interpolated, substitutes the $q_{i}$ for the $c_{i}$, and runs the row reduction step (see Section 2.2), that is, computes

[TABLE]

The transition then stores the actual result in the file-system and produces an output token which contains both a reference to the result and the point $q$. The stored data will be used later in the normalization step of the interpolation (see above).

Transition degree:

This generates a substitution point $p^{(j)}\in\mathbb{Z}^{\{1,\dots,j-1,j+1,\dots,k\}}$ for each $1\leq j\leq k$ yielding a matrix

[TABLE]

over the field $\mathbb{Q}(c_{j})$. After applying the row reduction, $M^{(j)}_{\text{red}}$ can be used to determine degree bounds for the numerator and denominator of each entry of the final result $M_{\text{red}}$ as a polynomial in $c_{j}$.

For $j\leq r$, we need a global degree bound to determine the number of interpolation points. We thus take the maximum of all numerator and denominator degrees of entries of $M^{(j)}_{\text{red}}$, and store these as a vector in $\mathbb{N}_{0}^{\{1,\dots,r\}}$, which is put on the place $d_{v}$.

If $j>r$, two integer matrices will be produced, which store the degrees of the numerators and denominators of each entry of $M_{\text{red}}$, respectively. This information will be used later to filter out bad interpolation points, that is, points at which polynomial cancellation occurs (see Appendix A). The result is stored in the file system and a token with a reference to the result is put on the place $d_{m}$.

Note that degree is in fact modeled by a sub-Petri net which behaves in a hierarchical manner as a transition. In practice, we actually compute multiple matrices $M^{(j)}$ per parameter to reduce the probability of a bad point producing wrong degree bounds.

Transition points:

This transition takes the degree data in $d_{v}$ and initializes the interpolation tree described in Section 4.2 and depicted in Figure 4. This, in turn, produces the corresponding set of interpolation points, which are put as separate tokens on the place $p$.

Transition reduce:

This transition consumes a point $p^{\prime}\in\mathbb{Z}^{\{1,\dots,r\}}$ from the place $p$ and computes

[TABLE]

The resulting matrix together with its interpolation point are put on the place $m$. Since reduce performs parameter substitutions in rational function expressions, the computation may fail due to division by zero. If this happens, $m$.valid is set to false, otherwise it is set to true.

Transition replace failure:

An input token for which $m$.valid is false is consumed by the transition replace failure, which marks the respective interpolation point as failed in the interpolation tree. If necessary, the interpolation tree is extended by additional interpolation points, which are also put on the place $p$.

Transition normalize:

An input token for which $m$.valid is true is consumed by the transition normalize. This transition reads $M_{\text{ref}}$ and multiplies the input matrix referenced by $m$ with a suitable constant factor. It also compares the entries with the degree matrices in $d_{m}$ to identify bad interpolation points. The result is put on the place $n$. If the corresponding point was bad, n.valid is set to false, otherwise to true.

Transition replace invalid:

For an input token for which n.valid is false, the transition generates new interpolation points in a fashion similar to that in replace failure.

Transition store normalized:

For an input token for which n.valid is true, the transition marks the corresponding interpolation point as successful in the external storage. If enough interpolation points for a given parameter have been marked as successful, the storage produces a token on place $i$, which triggers the respective interpolation. If the point $(p^{\prime}_{1},\dots,p^{\prime}_{r})$ triggers the interpolation (which will then use further points of the form $(p^{\prime}_{1},\dots,p^{\prime}_{r-1},p^{\prime\prime}_{r})$), the result of the interpolation will be associated to the point $(p^{\prime}_{1},\dots,p^{\prime}_{r-1})$ in the interpolation tree. If there are not yet enough interpolation points, the transition produces a token which only contains $i$.valid with value false.

Transition discard:

This transition discards tokens with $i$.valid equal to false.

Transition interpolate:

Tokens with $i$.valid equal to true are consumed by this transition, which then retrieves the references to the input data for the interpolation from the interpolation tree, loads the respective data from the file system, and executes the interpolation. If (in the above notation) the token holds $(p^{\prime}_{1},\dots,p^{\prime}_{r-1})$, then for $(d_{v})_{r}+1$ many points the corresponding row reduced matrices are retrieved from the storage. Note that due to the tree structure of the interpolation tree, all these points must have the first $r-1$ coordinates equal to $(p^{\prime}_{1},\dots,p^{\prime}_{r-1})$. The interpolation is then performed entry-wise as explained in Appendix A.

Transition store interpolated:

This transition marks the current point $(p^{\prime}_{1},\dots,p^{\prime}_{r-1})$ in the interpolation tree as processed. If $r>1$, just like in store normalized, the transition produces an interpolation token for the next parameter. If $r=1$, we have arrived at the final result, and a token with $i$.valid equal to false is produced, which will then be discarded.

The Petri net contains additional infrastructure (not described here) which terminates the execution once no tokens exist any more on the places $i$ and $p$.

4.4 Parallel timings

To illustrate the efficiency of our approach, we consider the cut $\{1,3,4,5\}$ of the double pentagon diagram (see Section 6 for a discussion of all possible cuts). Choosing this particular cut, which is less complex than others, our computations finish even when only a small number of cores is involed. This is necessary to analyze the scaling of our algorithm. In Table 1, we give timings for different numbers of cores. All timings are in seconds, taken on the high performance compute cluster at the Fraunhofer Institute for Industrial Mathematics (ITWM). Each compute node provides two Xeon E5-2670 processors, which amounts to $16$ cores333Hyperthreading is disabled. running at a base clock speed of 2.6 GHz. Each node has 64 GB of memory. For all runs with more than 15 cores, on each node we ran 15 compute jobs and one job for interfacing with the storage system. Since the storage jobs use negligible computation time, we omit them from the CPU core count when determining speedup and efficiency.

Apart from the running time $T(n)$ of the algorithm on a total of $n$ cores, we also give the speedup $S(n)=\smash{\frac{T(1)}{T(n)}}$ and the efficiency $E(n)=\smash{\frac{T(1)}{nT(n)}}$, which measure how “well” the algorithm parallelizes with increasing core counts. Note that the single-core timing is somewhat special: As experiments have shown, the performance per core decreases with the number of cores used on a given node. This effect has been investigated in BDFPRR (see in particular (BDFPRR, , Figure 5)). Thus, for the analysis of the expected run-time below, we rather consider the relative speedup and efficiency with respect to the 15-core timing. This in particular makes the assumption that the 15-core speedup is 15.

The saw-tooth shape of the efficiency graph in Figure 6 (and the corresponding behavior in the timing and speedup graphs) is due to the fact that the number of reductions to execute is usually not divisible by the number of cores utilized. Since in our test problem approximately $450$ reductions are required to enable the final interpolation, the running time of the full algorithm is roughly

[TABLE]

This effect can be avoided by a more fine-grained structuring of the problem (for instance by interpolating more parameters). Note, however, that increasing the number of processes in this way will lead to more overhead via inter-process communication and disk accesses. Thus, dividing the algorithm into very small parts may in fact slow down the overall computation.

Figure 6 also depicts the ideal expected runtime, speedup and efficiency. These ideal graphs stem from the simple assumption, called Amdahl’s law, that an algorithm can be divided up into a part that is ideally parallelizable and a part which is not parallelizable at all. Denoting the parallelizable fraction by $f$, the expected runtime $T_{\text{ideal}}(n)$ on $n$ cores is not $\frac{T(1)}{n}$, but rather

[TABLE]

which yields the ideal speedup and efficiency

[TABLE]

Using the experimental values for 15 and 30 cores, we arrive at a value $f\approx 0.999748$, that is, only 0.025 % of the algorithm is not parallelizable.

As we can see, the ideal curves give a fairly tight bound on the actual timings, at least in the cases where the core count is properly aligned to the number of reductions. This indicates that our approach for parallelization not only provides an automatic and fast solution to a tedious and complicated task, but stays highly efficient even when used with a large amount of computing power.

5 IBP conversion between different integral bases

It is well known that the IBP coefficient size may vary significantly if we choose different master integral bases. We prefer the IBP reduction to a uniformly transcendental (UT) basis as introduced in Henn:2013pwa , for several reasons: a) The differential equations satisfied by a UT integral basis have a particularly simple form Henn:2013pwa which allows for the integrals to be solved analytically in terms of polylogarithms. There is also evidence that for numerical computations, a UT basis is more convenient to evaluate444Private communication with Yang Gang.. So the IBP reduction to a UT basis greatly simplifies the amplitude computations after the IBP reduction. b) We observe that, in the case of the double pentagon, the IBP coefficients in a UT basis are significantly simpler than those in a traditional Laporta basis. This makes the IBP relations easier to use.

In practice, we consider special forms of UT bases, the so-called dlog bases, which will be introduced in the next subsection.

5.1 Dlog bases and the dlog algorithm

We say that a Feynman integral is a dlog integral if its integrand

[TABLE]

with $R(x_{1},...,x_{n})$ a rational function in $x_{1},...,x_{n}$, can be expressed in dlog form Bern:2014kca , that is, it can be written as a linear combination of type

[TABLE]

with rational functions $f_{i,j}$ in $x_{1},...,x_{n}$. This is only possible if the integrand has at most simple poles, including points at infinity. For example, both forms $\frac{dx}{x^{2}}$ and $dx$ admit no dlog form because of the double poles at zero respectively infinity.

The coefficients $c_{i}$ in equation (5.2) are called leading singularities ArkaniHamed:2010gh . For Feynman integrals, that are not of the elliptic type, they are in general algebraic functions of the external variables. By choosing an appropriate parametrization of the external variables, the leading singularities are typically rational functions. This is, in particular, true for the two-loop five-point integrals that are discussed in the next section. The leading singularities can also be understood as integrals over the original integrand where the integration contour is localized around the poles of the integrand. Leading singularities and the integrals integrated on the real contour have analytic properties in common. So, integrals with leading singularities that are just constant numbers are particularly useful, most importantly because they fulfill differential equations in the canonical form Henn:2013pwa . This implies that they have the property of uniform transcendental weight, which means that if the series is expanded in $\epsilon$, the parameter of dimensional regularization, the coefficients have homogeneous transcendental weight and the weight increases by one for each order in $\epsilon$.

Next, we recall from WasserMSc how to transform a given integrand into dlog form, in case this is possible. Given an integrand in $n$ integration variables, we choose, if possible, one variable $x$ that is linear in all denominator factors and do a partial fraction decomposition while treating all other variables as constants. In this way, we obtain a sum of integrands of the form

[TABLE]

where $\Omega$ is an $(n-1)$-form, independent of $x$, and $a$ is a polynomial that may depend on the other integration variables. Then we iterate this procedure taking $\Omega$ as our new integrand until no integration variables are left. If in any intermediate step a pole of degree two or higher is encountered, then the integrand does not admit a dlog form. There are cases where no variable exists that is linear in all denominator factors. One way to proceed in such a case is to make a variable transformation such that at least one of the new variables is linear in all denominator factors.

The algorithmic approach of this section was used in Chicherin:2018mue and Chicherin:2018old to construct a complete basis of dlog master integrals with constant leading singularities for all two-loop five point integral families. The denominator structure for each integral family is given by the propagators. To construct the dlog integrals we make a general numerator ansatz. We write the numerator as a linear combination of terms that are products of inverse propagators and irreducible scalar products. Each term is multiplied by a free parameter, and by applying the algorithm to this general integrand, we can determine values of the free parameters such that the integrand has a dlog form and constant leading singularities. In this way, we obtain a set of dlog integrals that form a basis of dlog master integrals.

In general, the dlog algorithm can be applied only in a dimension that is an integer number, which we choose to be four. The loop momenta are very conveniently parametrized using spinor helicity variables as in Bern:2014kca . Although this parametrization can be very useful, it also has its limitations as soon as the numerator has terms that vanish in dimension four, but which are non-zero in generic dimension $D$. In such cases, an extended approach as in Chicherin:2018old using the Baikov parametrization can be applied.

5.2 IBP reduction with a dlog basis

Given a dlog basis, we discuss the IBP reduction in two settings:

When both the IBP coefficients in the Laporta basis and the dlog basis are needed, we first compute the reduction in the Laporta basis $I$ with our module intersection and GPI-Space reduction algorithm,

[TABLE]

where $F$ is the list of target integrals as a column vector. Then we reduce the dlog basis $\tilde{I}$ to the Laporta basis $I$,

[TABLE]

Note that since the dlog basis construction has a restriction on the numerator degree, this reduction is usually easy. Terms exceeding the allowed numerator degree have double poles at infinity. This can be seen by inverting the loop momenta $k_{i}^{\mu}\rightarrow k_{i}^{\mu}/k_{i}^{2}$. Using our Singular RREF code, with a good pivot strategy, we can analytically find the inverse $T^{-1}$. The matrix product $AT^{-1}$ contains the coefficients of an IBP reduction to the dlog basis.

We remark that the product $AT^{-1}$ can be difficult to calculate even if $T^{-1}$ has a relative small size. Instead of computing the product directly, we again use the semi-numerical approach, setting several of the kinematic values to be integers, computing the product several times, and then using our interpolation program to get the fully analytical matrix product $AT^{-1}$. This is again implemented using our Singular-GPI-Space framework. 2. 2.

When only the IBP coefficients in a dlog basis are needed, we apply our semi-numerical reduction method to a set of numeric IBP coefficients in the Laporta basis. Instead of interpolating these coefficients, we use the semi-numeric points to interpolate the product $AT^{-1}$, not calculating the analytic form of $A$.

In the next section, we illustrate our approach by considering a non-trivial example, the two-loop five-point nonplanar double pentagon diagram. This includes the IBP generation via the module intersection method, the massively parallel reduction of the IBP system and the basis conversion.

6 The two-loop five-point nonplanar double pentagon example

In this section, we illustrate our IBP reduction method by applying it to a nontrivial example, the two-loop five-point nonplanar double pentagon. Note that a symbolic UT basis for this example was found in Abreu:2018aqd ; Chicherin:2018yne . Furthermore, UT bases in terms of polylogarithm functions for the double pentagon and other two-loop five-point nonplanar massless integral families were analytically calculated in Chicherin:2018old .

For the diagram in Figure 7, we chose the following labeling for the propagators:

[TABLE]

where the $l_{i}$ represent the loop momenta, the $p_{i}$ represent external momenta, and $p_{i\cdots j}=\sum_{i}^{j}p_{i}$. The first $8$ propagators represent the topology and the last three ones the irreducible scalar products.

This is a complicated integral family for IBP reduction, due to the number of independent scalars, which are s12, s23, s34, s45, s45, s15 and the spacetime dimension D, and due to the nonplanar topology with two pentagons inside. We demonstrate our method by reducing the $26$ integrals listed in Figure 8

to a master integral basis in the fashion of Laporta. Furthermore, we convert the IBP coefficients to the coefficients of a dlog basis given in Chicherin:2018old . In this base change, we observe a significant coefficient size reduction.

6.1 Module intersection with cuts

First, we use Azurite Georgoudis:2016wff to find an integral basis. Without considering symmetries, there are $113$ irreducible integrals, and with symmetries, there are $108$ master integrals. Note that due to the number of master integrals, this IBP reduction is significantly more complicated than the reduction of the hexagon-box diagram in Boehm:2018fpv , which has only $73$ master integrals.

We can then construct the set of spanning cuts of this integral family. Each spanning cut corresponds to a “non-collapsible” master integral Georgoudis:2016wff . There are $11$ spanning cuts (without considering symmetries),

[TABLE]

where the numbers indicate the propagators on cut. For example, $\{3,4,7\}$ means that $D_{3}\mapsto 0,D_{4}\mapsto 0,D_{7}\mapsto 0$.

For each cut, we can apply our module intersection method to generate IBPs without double propagators. In Boehm:2018fpv , the IBPs are generated only from the top sector. Here, for each cut, we generate the IBPs from both the top sector and lower sectors. For example, for the cut $\{1,5,7\}$, we consider the $32$ sectors supported on cut $\{1,5,7\}$, and compute all the module intersections. This approach will generate more IBPs but make the IBP system more block-triangular and easier for linear reduction.

With the degbound option in Singular, it is easy to generate all the module intersections. For this integral family, choosing the degree bound $5$, and using one core for each cut, it takes less than $5$ minutes in total to solve all the module intersection problems analytically. Later on, by finite-field methods, we find that with this choice of degree bound, we obtain sufficiently many IBPs for our problem.

After generating the IBPs, we use the two-step trimming process described in Section 2.2 to select necessary IBPs for our targets. This computation is via finite-field methods and powered by the package SpaSM.

Note that different cuts can support the same master integral. We then speak of a cut overlap. For example, the integral $I[1,1,1,1,1,0,1,1,0,0,0]$ is supported by both cuts $\{1,5,7\}$ and $\{2,4,8\}$. The IBP reductions on these two cuts should give the same coefficients. To avoid redundant computations, we can apply the master integral removal method Chawdhry:2018awn , setting $I[1,1,1,1,1,0,1,1,0,0,0]\mapsto 0$ either in cut $\{1,5,7\}$ or in cut $\{2,4,8\}$. Clearly there are many different removal strategies for master integral overlapping in different cuts, and different strategies result in different performances. In our computational problem, we find that the cuts $\{1,6,8\}$ and $\{2,4,8\}$ are relatively “difficult” cuts for IBP reduction. Hence, we set as many master integrals as possible in $\{1,6,8\}$ and $\{2,4,8\}$ to zero, and later on recover the remaining master integral coefficients from other cuts via cut overlap.

We compute the module intersections analytically. For the purpose of linear reduction, we further set

[TABLE]

to dehomogenize the IBP relations and speed up the computation. The $s_{12}$ dependence can be recovered in the final step.

The resulting IBPs are summarized in Table 2. Note that for the cut $\{1,6,8\}$, there are $1203$ independent relations and $1205$ integrals after applying the idea of Chawdhry:2018awn to set most master integrals supported on the cut $\{1,6,8\}$ to zero. As a result we only have to compute just two master integral coefficients.

6.2 IBP reduction

We apply our reduction method via Singular and GPI-Space to reduce the linear systems in Table 2. We use a semi-numeric approach, choosing $c_{4}$, $c_{5}$ and the space-time dimension $D$ to be symbolic, and compute the linear reduction with integer-valued $c_{2}$ and $c_{3}$.

By a linear reduction with $c_{2}$ (respectively $c_{3}$) symbolic and all the other parameters numeric, we easily determine the maximal degree of $c_{2}$ (respectively $c_{3}$) in the reduced IBP relations. The degrees are listed in Table 2 as $d_{2}$ and $d_{3}$, respectively. From this information, we get the minimal number $(d_{2}+1)\times(d_{3}+1)$ of semi-numeric computations for interpolating the analytic reduction result. For example, for the cut $\{1,5,7\}$, we need to run semi-numeric computations at least $506$ times. Of course, the cuts exhibit different running times when performing the reductions: For instance, cut $\{1,3,4,5\}$, which we already considered as an example in Section 4.4, is the easiest in terms of running time, taking only about $11$ minutes when using $384$ CPU cores. In contrast, the cut $\{3,4,8\}$ is much more complex: its reduction took $12$ hours and $21$ minutes, using $384$ cores.

After getting the analytic reduction of all the cuts, we merge them to get the full IBP reduction to a $113$-integral basis. Furthermore, we apply the symmetries, obtained from Georgoudis:2016wff ,

[TABLE]

to reduce the $26$ target integrals to a $108$-integral Laporta basis $I$. We note that the resulting file is large, with a size of $\sim 2.0$ GB on disk.

By setting all Mandelstam variables to integers, we have verified that our result is consistent with FIRE6 Smirnov:2019qkx .

6.3 IBP coefficient conversion to a dlog basis

In this subsection, we discuss converting the IBP coefficients for the Laporta basis to the IBP coefficients of the dlog basis found in Chicherin:2018old .

For this conversion, we again use the semi-numeric approach, taking integer-valued $c_{2}$, $c_{3}$, and symbolic $c_{4}$, $c_{5}$ and $D$, converting the coefficients and then interpolating. It is easy to determine that the coefficients in the dlog basis have the following maximal degrees for $c_{2}$ and $c_{3}$, respectively,

[TABLE]

By comparing with Table 2, where $d_{2}$ can be as high as $35$, we find that the maximal degree drops. For the basis conversion, we carry out a semi-numeric matrix multiplication with subsequent interpolation using Singular and GPI-Space.

After the computation, we see that the IBP reduction coefficients of Figure 8 in this dlog basis have size $480$ MB on disk, which shows a significant $76\%$ reduction of the IBP coefficient size compared to what we have for the Laporta basis. On the other hand, if only the IBP reduction coefficients in the dlog basis are needed, we can skip the interpolation for the Laporta basis IBP coefficients, and directly convert the intermediate numerical results to dlog basis IBP coefficients. Because of the maximal degree drop, this shortcut reduces the required number of semi-numeric computations.

For convenience, we also provide the IBP coefficients in the dlog basis, with the $s_{12}$ scalar recovered. All these analytic results can be obtained via the links presented in the introduction of this paper. Note that all files provided under the links contain $26\times 108$ matrices. For each matrix, the entry in the $i$th row and $j$th column is the corresponding IBP coefficient for the $i$th target integral in Figure 8, expanded on the $j$th master integral. The Laporta basis and the dlog basis are included in the auxiliary files of this paper.

7 Summary

In this paper, we present our powerful new IBP reduction method, which is based on computational algebraic geometry powered by the computer algebra system Singular in conjunction with the taskflow management system GPI-Space. Our method is suitable for large scale IBP reduction problems with complicated Feynman diagrams and multiple variables. We demonstrate the power of the new method by the analytic two-loop five-point nonplanar double pentagon IBP computation. The computational result has been cross-checked numerically using state-of-the-art IBP programs.

Our method is flexible and can be adapted in various different scenarios:

Modern methods for amplitude computation often follow the approach of numerically or semi-numerically calculating the IBP relations in order to interpolate the amplitude coefficient under consideration directly, instead of interpolating the analytic IBP relations. Our method can efficiently compute the reduced numeric or semi-numeric IBP relations and, hence, perfectly fits into this purpose. 2. 2.

Our module intersection method can also be used for integrals with double propagators or multiple-power propagators since this IBP generating method avoids the increase of propagator exponents and significantly reduces the size of the IBP system. 3. 3.

Although our method is currently based on semi-numerical parallelizations with integer-valued numerics, it clearly can be extended to finite-field linear reduction, if necessary. 4. 4.

More generally, our linear reduction parallelization method can be used for computational problems other than IBP reduction. For example, in recent years, it was found that the Bethe Ansatz equation of integrable spin chains can be analytically computed by algebraic geometry methods Jacobsen:2018pjt ; Jiang:2017phk . Often, this involves large-scale linear algebra computations with symbolic parameters, and our parallelization via the Singular-GPI-Space framework can greatly speed up the computation. We also expect that our reduction method can be used more generally for Gröbner basis computations with parameters.

In the future, we will develop our code into an automated software package, powered by Singular and GPI-Space, for solving large-scale IBP or amplitude problems. The possible simplification of IBP coefficients in a UT/dlog basis will be further investigated. We expect that this method will provide a boost for the current NNLO precision computations.

Acknowledgments

YZ thanks Johann Usovitsch for many enlightening discussions on IBP reduction algorithms and their implementation, and Johannes Henn for cutting-edge techniques for finding uniformly transcendental and dlog integrals. We acknowledge Kasper Larsen for his work in the starting stage of the analysis of the two-loop nonplanar five-point double pentagon diagram. We thank Simon Badger, Rutger Boels, Christian Bogner, Hjalte Frellesvig, Gudrun Heinrich, Harald Ita, Stephan Jahn, David Kosower, Roman Lee, Hui Luo, Yanqing Ma, Andreas von Manteuffel, Pierpaolo Mastrolia, Alexander Mitov, Erik Panzer, Tiziano Peraro, Robert Schabinger, Hans Schönemann, Alexander Smirnov, Vladimir Smirnov, Peter Uwer, Gang Yang and Mao Zeng for their help on this research project.

The research leading to these results has received funding from the Swiss National Science Foundation (Ambizione grant PZ00P2 161341), from the National Science Foundation under Grant No. NSF PHY17-48958, from the NSF of China with Grant No. 11235010, from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 725110). The work of AG is supported by the Knut and Alice Wallenberg Foundation under grant #2015-0083. The work of DB, JB, and WD was supported by Project II.5 of SFB-TRR 195 Symbolic Tools in Mathematics and their Application of the German Research Foundation (DFG). The authors would also like to express a special thanks to the Mainz Institute for Theoretical Physics (MITP) of the Cluster of Excellence PRISMA+ (Project ID 39083149) for its hospitality and support during the “Mathematics of Linear Relations between Feynman Integrals” workshop.

Appendix A Rational function interpolation

In this appendix, we introduce our simple approach to rational function interpolation. Although this algorithm is rather straight-forward compared to other more involved techniques available (see, for example, zippel ; bentiwari ; KLW ), we have found that it is more suitable for our setup. The idea is to heuristically convert a rational function interpolation problem to a polynomial function interpolation problem.

We focus on a general computational process with symbolic variables $x_{1},\ldots x_{k}$ which would give the final result as a rational function

[TABLE]

Here, $F$ and $G$ are integer valued polynomials with $\gcd(F,G)=1$. Suppose we are in a situation where it is difficult to perform this process with all parameters $x_{1},\ldots,x_{k}$ symbolic, but where it is feasible to obtain results when some of the parameters are set to be general integer values, and only the remaining parameters are symbolic. Then our approach is to repeat such semi-numeric computations sufficiently often, and interpolate to get (A.1).

To be specific, in such a situation, and for a fixed $k_{1}<k$, we refer to the computation with $x_{k_{1}+1},\ldots,x_{k}$ symbolic, and the other parameters substituted by random integers,

[TABLE]

as the $i$th semi-numeric computation. We write the result of the $i$th semi-numeric computation as

[TABLE]

Note that, although computer algebra software can cancel the fraction to get polynomials $f_{i}$ and $g_{i}$ with $\gcd(f_{i},g_{i})=1$, the relation between $G$ and $g_{i}$ is not clear a priori. For example, it may happen that

[TABLE]

Similarly, we may have that

[TABLE]

The reason for this is that after taking integer values $a_{1}^{(i)},\ldots,a_{k}^{(i)}$ for the first $k$ variables, there can be additional cancellations between $F$ and $G$. This phenomenon makes the direct polynomial interpolation of the $g_{i}$ and $f_{i}$ inapplicable.

We solve this cancellation problem in a heuristic way.

First, we compute a “reference” result with symbolic $x_{1},\ldots,x_{k}$ and random integer values for the other parameters, $x_{k_{1}+1}\mapsto b_{k_{1}+1},\ldots,x_{k}\mapsto b_{k}$:

[TABLE]

where $p$ and $q$ are integer valued polynomials with $\gcd(p,q)=1$. Generally, except for a statistically very small set of points $b=(b_{k_{1}+1},\ldots,b_{k})$, we can assume that the two polynomials $F(x_{1},\ldots,x_{k_{1}},b_{k_{1}+1},\ldots,b_{k})$ and $G(x_{1},\ldots,x_{k_{1}},b_{k_{1}+1},\ldots,b_{k})$ are coprime. Then, since

[TABLE]

by the unique factorization domain (UFD) property,

[TABLE]

for some integer $\mu$. Note that our algorithm relies only on the existence of such a $\mu$. 2. 2.

Second, we do a “majority vote” selection from among the semi-numeric results: For each of the polynomials, we record the leading exponent in $x_{k+1},\ldots,x_{k}$, and determine the most frequently occurring exponent $\mathbf{r}$. Then we drop all semi-numeric result whose exponent is not equal to $\mathbf{r}$. This step ensures that the gcd of $F(a_{1}^{(i)},\ldots,a_{k}^{(i)},c_{k_{1}+1},\ldots,c_{k})$ and $G(a_{1}^{(i)},\ldots,a_{k}^{(i)},x_{k_{1}+1},\ldots,x_{k})$ is just an integer instead of a non-constant polynomial in $x_{k_{1}+1},\ldots,x_{k}$. Again by the UFD property,

[TABLE]

where each $\lambda_{i}$ is an integer. Setting $x_{1}\mapsto a_{1}^{(i)},\ldots,x_{k_{1}+1}\mapsto a_{k_{1}+1}^{(i)}$ in (A.8) and $x_{k_{1}+1}\mapsto b_{k_{1}+1},\ldots,x_{k}\mapsto b_{k}$ in (A.10), we determine that

[TABLE]

3. Define

[TABLE]

Then interpolate both polynomials $\tilde{f}_{i}(x_{k_{1}+1},\ldots,x_{k})$ and $\tilde{g}_{i}(x_{k_{1}+1},\ldots,x_{k})$ for the selected semi-numeric points by standard polynomial interpolation algorithms, say by using Newton polynomials. From (A.12), we see that the resulting two polynomials are

[TABLE]

so that the ratio of the polynomials gives the desired fraction $F/G$, while the integer factor $\mu$ is canceled out.

So with the compensation factors $\lambda_{i}/\mu$, we can use simple polynomial interpolation algorithms to interpolate rational functions at the extra cost of computing only one additional reference point. If the degree of $F$ and $G$ in $x_{j}$ in is $d_{j}$ ($1\leq j\leq k_{1}$), for a general choice of the interpolation points, we need to compute the semi-numeric result

[TABLE]

times, plus one computation for the reference point.

The resulting algorithm is implemented using Singular in conjunction with GPI-Space.

In practice, this algorithm is extendable in many ways. For example, instead of splitting the variables $x_{1},\ldots x_{k}$ into two groups $(x_{1},\ldots,x_{k_{1}})$ and $(x_{k_{1}+1},\ldots,x_{k})$, we can split the variables in more groups and use our algorithm recursively. The algorithm can also be combined with finite field reconstruction.