I'm using VEGAS integration, specifically the GSL implementation, for some QCD calculations, and I've been investigating the behavior of the algorithm for various numbers of iterations in an attempt to get more accurate results. The value produced by the integration routine changes rather drastically with the number of iterations, as shown in this plot:

plot of value increasing with iterations

The results shown are for one specific random seed, but they are indicative of other results obtained by running many instances of the program with different random seeds and averaging; specifically, the value with a million iterations is significantly larger than that with a hundred thousand iterations, as shown in this plot (which includes the above value, plus a few other terms whose values are nearly constant):

box plot of statistical results

My best guess as to the cause is that there are peaks in the function being integrated which have a characteristic size around $V/10^6$, with $V$ being the integration volume (in this case, an eight-dimensional hypercube). These peaks would then be detected when sampling with a million points, but sampling with only a hundred thousand points would have a good chance of missing them. That would explain why some of the values at $10^5$ iterations are larger than the values with $10^6$ iterations.

Preliminary question: am I on the right track with this explanation?

Assuming I am, that brings up the question of how to know when I can safely stop increasing the number of iterations. For example, is there some sort of analysis I can apply to the analytic form of the function to pick out the scale of features which will significantly affect the result of the integration?

The integral in question, in case it helps:

$$\begin{multline} \int_{\tau}^1\mathrm{d}z\int_{\tau/z}^1\mathrm{d}\xi\iiint\mathrm{d}^2\vec{q}_1\mathrm{d}^2\vec{q}_2\mathrm{d}^2\vec{q}_3 \\ \frac{1}{1 - \xi}\biggl[\frac{\frac{\tau}{z\xi}g\bigl(\frac{\tau}{z\xi}\bigr)D(z)}{z^2}\frac{(1 - \xi(1-\xi))^2}{\xi^2}F\biggl(\frac{\vec{k}}{\xi} + \vec{q}_3 - \vec{q}_2, z\biggr) - \frac{\frac{\tau}{z}g\bigl(\frac{\tau}{z}\bigr)D(z)}{z^2}F\bigl(\vec{k} + \vec{q}_3 - \vec{q}_2, z\bigr)\biggr] \\ \times F\bigl(\vec{k} + \vec{q}_3 - \vec{q}_1, z\bigr)F(\vec{q}_3, z)\frac{\vec{q}_1\cdot\vec{q}_2}{q_1^2q_2^2} \end{multline}$$ where $F(\vec{q}, z) = \frac{z^{0.288}}{C\pi}\exp\bigl(-q^2 z^{0.288}/C\bigr)$ with $C$ a constant.


The integral in question has a singularity at $\xi=1$, and this singularity is at the border of the integration domain. You have to handle singularities analytically or remove them by some other mean, before you can use numerical integration, and VEGAS integration is no exception here.

You probably have some idea why the value of the integral is well defined and finite despite this singularity. Such an idea normally also suggests a method how the singularity could be removed or handled analytically.

  • $\begingroup$ Actually the singularity has already been removed. That's the purpose of the second term in square brackets. $\endgroup$ – David Z Apr 13 '15 at 4:55
  • $\begingroup$ To whoever downvoted this wrong answer: This was exactly the right thing to do. I won't delete this answer, because I'm not a friend of hiding failed solution attempts. However, it is still wrong, and hence should be downvoted. $\endgroup$ – Thomas Klimpel Aug 11 '15 at 6:49

Your Answer

