Vectorial integration

[alecandido]

The problem

In order to compute the singlet operators eko needs four elements: S_qq, S_qg, S_gq, S_gg (in get_ker(k, l) only the k, l element is provided).

These elements will depend on the anomalous dimensions, but in order to compute each of them all the 4 anomalous dimensions should be evaluated, cf. get_gamma_singlet_0, called by get_Eigensystem_gamma_singlet_0.

But since the eigenvalue decomposition has to be performed N by N (Mellin) is not possible to integrate before the 4 anomalous dimensions and then rotate to the operators, but the single evaluation for a specific N will involve the decomposition, and so the calculation of all of the 4 anomalous dimensions.

What is happening so is that for a given N the anomalous dimensions are evaluated every time for each of the 4 integration of S_qq, S_qg, S_gq, S_gg, and in the integration process evaluating the integrand will take most of the time (almost the total amount, we believe).

A significant improvement would be integrate simultaneously the four elements, indeed when the singlet-gluon sector is relevant (in fact: always) all of the four elements are needed, and in this way all of the four elements produced by get_Eigensystem_gamma_singlet_0 are used, instead of taking one and discarding 3 every time.

Considering also the time required for the non-singlet calculation this will lead to a speed up naïvely estimated by a factor of 3.4=(4*4+1)/(4+1), independently of the complexity of the kernel (LO, NLO, NNLO, ... with whatever method used).

Naïve meaning
If not any other error is present at least this estimate is rough mainly because of the different complexity of the integration of the different elements.
Indeed if any relevant cancellation will smooth some of the elements the convergence will be much faster than for the other elements, because less evaluation are needed.
However in the worst case there is no gain and no loss in performance, since the time is spent mostly in integrating the most complex element, but all the evaluation done are already done by the current implementation, and in this way the other 3 entries will be used to improve the precision of the simplest elements.

Useless solution

scipy.integrate subpackage has a function called scipy.integrate.quad_vec, but it is currently implemented in python (while scipy.integrate.quad currently used is a binding for QUADPACK fortran implementation).

Since the latter is compiled switching to quad_vec it is likely to have no performance benefit.

Proposed solutions

• drop scipy.integrate.quad and fall back to a sum
  • this will reduce the speed up, but it is quick to implement and benchmark
• implement scipy.integrate.quad_vec in a C/Fortran compiled module of eko
  • this will get the theoretical speed up, but we will increase the complexity of package management (so we will loose in development time)
• implement scipy.integrate.quad_vec in a numba compiled function, inside eko
  • in principle can obtain the same speed of the former, but it will be easier to manage and distribute
• implement scipy.integrate.quad_vec in C/Fortran and pull request into scipy
  • slower development for the function itself (not familiar environment, external people slighlty involved), optimal solution for not polluting eko and keeping the delegation to an external utility

[scarlehoff]

At some point I talked with @scarrazza about giving this to a student as it is not something extremely complicated but can be nice as a first experience with numerical integrations. And would be useful beyond eko because as you say the scipy implementation is really bad (performance wise).

There is also an extra thing we can do, related with issue #56. One can create a Cuda-quad_vec (and write the cuda-kernels necessary to make the whole of the integrand into a kernel). I don't know how much it will improve (can't be known until we try!)
Of course, this is 100% a side project, but I thought it was worth to mention.

[alecandido]

Thank you very much @scarlehoff, I don't have much experience with Cuda nor with GPU in general, but with generic parallelization yes, at least a bit, and if possible it would be great.

Since as you said it's quite a separate project (even if will be completely useful for eko) I would suggest that if at any time you have some man-power it's fine to make it (it's even possible that we will have some more man-power...).

We will take care for sure of the integration part, but until some external help won't be available I don't think we will have time enough to spend on this, unfortunately. Just because at the time I wrote this NLO was running in 10-30 minutes, but now it's already running in less than 1, at least often, so it's definitely not a priority anymore (but maybe it will be very useful for NNLO).

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

[alecandido]

I believe this is going to be critical for QED, since the coupled sector is going to increase. I would assign more priority.

[felixhekhorn]

I'm not sure - there are many more tasks with higher priority ... (like parallelize)

[felixhekhorn]

@alecandido another candidate for https://github.com/N3PDF/eko/wiki/List-of-Possible-Future-Features ?

[alecandido]

Approved