# Multi-dimensional integral fails to converge using GSLMC VEGAS

Hello,

I am using the `ROOT::Math::GSLMCIntegrator` with `kVEGAS` to integrate a multi-dimensional function. The function is a polynomial with negative exponentials, defined as ff6 in the code. The integration of ff6 gives the matrix elements of row number 6. This is a 10x10 correction matrix. I would like to calculate the integrals with an absolute tolerance of at least 1e-4 for all elements of the matrix. However, not even for an absolute tolerance of 1e-1 the integration converges. ROOT keeps running it forever.

The same integration if performed using `ROOT::Math::AdaptiveIntegratorMultiDim` with absolute tolerance of 1e-3 takes about 4 hours to complete. When I calculate the elements in rows 8, 9 and 10, each element takes around 10 hours with an absolute tolerance of 1e-3. With an absolute tolerance of 1e-4 it would take days to run a few elements. I tried substituting the exponentials by a series expansion, but the series only converge after 20 terms which would complicate even more the function being integrated. I did this because I have the impression that the exponential terms are the ones responsible for slowing down the integration.
Do you have an idea of why the MC integration is not converging? Is there another method I could try?

Kind regards,
Jeremias

``````double f(double x) {

TF1 *tf1 = new TF1("tf1","(**2*)/(-)**2*(exp(-*x)*((-)*x-1)+exp(-*x))*+",0.0,100.0); // is lambda and  stands for beta.
tf1->SetParameters(0.154,0.0722,100.0,0.0001,0.0);
Double_t top=tf1->Eval(tf1->GetMaximumX())+100.0;
TF1 *tf2 = new TF1("tf2","tf1-",0.0,100.0);
TF1 *tf4 = new TF1("tf4","tf2/()",0.0,60.0);
tf4->SetParameter(0,tf2->Integral(0.0,60.0));

Double_t p0=tf1->GetParameter(0), p1=tf1->GetParameter(1), p2=tf1->GetParameter(2), p3=tf2->Integral(0.0,60.0);

delete tf1; delete tf2; delete tf4;

return ((pow(p0,2)*p1/pow(p1-p0,2))*(exp(-p0*x)*((p1-p0)*x-1.0)+exp(-p1*x))*p2)/p3;
}

//____________________________

double fint1(double t)
{
ROOT::Math::Functor1D wf(&f);
ig.SetFunction(wf);
Double_t tau=0.9;
Double_t fintval = ig.Integral(t,t+tau);

return fintval;
}
//_______________________________i=6-------------------------------------------------------

double ff6(const double *x, const double *par)
{
const int i=6;
Double_t tau=0.9;
Double_t T=60.0;
Double_t t2min=x+tau;
Double_t t2max=T-(i-2.0)*tau;
Double_t t2=(t2max-t2min)*x+t2min;
Double_t dt2=t2max-t2min;
Double_t t3min=t2+tau;
Double_t t3max=T-(i-3.0)*tau;
Double_t t3=(t3max-t3min)*x+t3min;
Double_t dt3=t3max-t3min;
Double_t t4min=t3+tau;
Double_t t4max=T-(i-4.0)*tau;
Double_t t4=(t4max-t4min)*x+t4min;
Double_t dt4=t4max-t4min;
Double_t t5min=t4+tau;
Double_t t5max=T-(i-5.0)*tau;
Double_t t5=(t5max-t5min)*x+t5min;
Double_t dt5=t5max-t5min;
Double_t t6min=t5+tau;
Double_t t6max=T-(i-6.0)*tau;
Double_t t6=(t6max-t6min)*x+t6min;
Double_t dt6=t6max-t6min;
Double_t ffval = (TMath::Factorial(par)/TMath::Factorial(par-i)) *f(x) *f(t2) *f(t3) *f(t4) *f(t5) *f(t6) *dt2 *dt3 *dt4 *dt5 *dt6 *pow(fint1(x) +fint1(t2) +fint1(t3) +fint1(t4) +fint1(t5) +fint1(t6),par-i);
return ffval;
}

//____________________________

double S6j(double j)
{
const int i=6;
Double_t tau=0.9;
Double_t T=60.0;
Double_t epsabs=1e-1;
ROOT::Math::WrappedParamFunction<> wf( &ff6, i, 1);
Double_t p = {j};
wf.SetParameters(p);
// Create the Integrator
ROOT::Math::GSLMCIntegrator ig(ROOT::Math::MCIntegration::kVEGAS);
// Set parameters of the integration
ig.SetFunction(wf);
ig.SetAbsTolerance(epsabs);
Double_t xmin[i] = {0.0,0.0,0.0,0.0,0.0,0.0};
Double_t xmax[i] = {T-(i-1.0)*tau,1.0,1.0,1.0,1.0,1.0};
Double_t Sjvalue =ig.Integral(xmin, xmax);
return Sjvalue;
}
void spencerIntegrals() {
ROOT::EnableImplicitMT(0);
cout << "Matrix element S66: " << S6j(6) << endl;
cout << "Matrix element S67: " << S6j(7) << endl;
}
``````

ROOT Version: v6-28-06
Platform: Red Hat Enterprise Linux 8.8; Memory 125.1 GiB; Intel Core i7-9800x CPU 3.8 GHz x16.
Compiler: C++ (GCC) 8.5.0 20210514 (Red Hat 8.5.0-18)

Hi Garcia,

In the integrand you define two functions. It is called many times and that can not be optimal. I put the definitions outside the integrand, however it is unhappy in the GSL call but I let you figure that out:

``````#include "Math/Functor.h"
#include "Math/Integrator.h"
#include "Math/WrappedMultiTF1.h"
#include "Math/WrappedParamFunction.h"
#include "Math/GSLMCIntegrator.h"
#include "TF1.h"

TF1 *tf1;
TF1 *tf2;

auto sqr(double x)
{ return x * x; }

double f_1(const double *x, const double *par) {
Double_t val = (sqr(par)*par)/sqr(par-par)*(exp(-par*x)*((par-par)*x-1)+exp(-par*x))*par+par;
return val;
}
double f_2(const double *x, const double *par) {
Double_t val = f_1(x,par)-par;
return val;
}

double f(double x) {

Double_t p0=tf1->GetParameter(0), p1=tf1->GetParameter(1), p2=tf1->GetParameter(2), p3=tf2->Integral(0.0,60.0);
return ((sqr(p0)*p1/sqr(p1-p0))*(exp(-p0*x)*((p1-p0)*x-1.0)+exp(-p1*x))*p2)/p3;
}

//____________________________

double fint1(double t)
{
ROOT::Math::Functor1D wf(&f);
ig.SetFunction(wf);
Double_t tau=0.9;
Double_t fintval = ig.Integral(t,t+tau);

return fintval;
}
//_______________________________i=6-------------------------------------------------------

double ff6(const double *x, const double *par)
{
const int i=6;
Double_t tau=0.9;
Double_t T=60.0;
Double_t t2min=x+tau;
Double_t t2max=T-(i-2.0)*tau;
Double_t t2=(t2max-t2min)*x+t2min;
Double_t dt2=t2max-t2min;
Double_t t3min=t2+tau;
Double_t t3max=T-(i-3.0)*tau;
Double_t t3=(t3max-t3min)*x+t3min;
Double_t dt3=t3max-t3min;
Double_t t4min=t3+tau;
Double_t t4max=T-(i-4.0)*tau;
Double_t t4=(t4max-t4min)*x+t4min;
Double_t dt4=t4max-t4min;
Double_t t5min=t4+tau;
Double_t t5max=T-(i-5.0)*tau;
Double_t t5=(t5max-t5min)*x+t5min;
Double_t dt5=t5max-t5min;
Double_t t6min=t5+tau;
Double_t t6max=T-(i-6.0)*tau;
Double_t t6=(t6max-t6min)*x+t6min;
Double_t dt6=t6max-t6min;
Double_t ffval = (TMath::Factorial(par)/TMath::Factorial(par-i)) *f(x) *f(t2) *f(t3) *f(t4) *f(t5) *f(t6) *dt2 *dt3 *dt4 *dt5 *dt6 *pow(fint1(x) +fint1(t2) +fint1(t3) +fint1(t4) +fint1(t5) +fint1(t6),par-i);
return ffval;
}

//____________________________

double S6j(double j)
{
const int i=6;
Double_t tau=0.9;
Double_t T=60.0;
Double_t epsabs=1e-1;
ROOT::Math::WrappedParamFunction<> wf( &ff6, i, 1);
Double_t p = {j};
wf.SetParameters(p);
// Create the Integrator
ROOT::Math::GSLMCIntegrator ig(ROOT::Math::MCIntegration::kVEGAS);
// Set parameters of the integration
ig.SetFunction(wf);
ig.SetAbsTolerance(epsabs);
Double_t xmin[i] = {0.0,0.0,0.0,0.0,0.0,0.0};
Double_t xmax[i] = {T-(i-1.0)*tau,1.0,1.0,1.0,1.0,1.0};
Double_t Sjvalue =ig.Integral(xmin, xmax);
return Sjvalue;
}

void spencerIntegrals() {

TF1 *tf1 = new TF1("tf1",f_1,0.0,100.0,4);
tf1->SetParameters(0.154,0.0722,100.0,0.0001);
Double_t top=tf1->Eval(tf1->GetMaximumX())+100.0;
TF1 *tf2 = new TF1("tf2",f_2,0.0,100.0,4);

//ROOT::EnableImplicitMT(0);
cout << "Matrix element S66: " << S6j(6) << endl;
cout << "Matrix element S67: " << S6j(7) << endl;
}
``````
1 Like

Hi Eddy,

Thank you for your reply. In my code I was actually fitting the data with `tf1` inside the `f` function. After applying your advice and fitting the function outside the integrand, my code got 10 times faster using `ROOT::Math::AdaptiveIntegratorMultiDim`. Now I am able to go down to a tolerance of 1e-5.
However, if I use the `ROOT::Math::GSLMCIntegrator` with `kVegas`, the integration still takes a very long time, I don’t know if it will eventually converge.

Best regards,
Jeremias

Just as a reference, this is a simplification of my code if anyone wants to try:

``````#include "Math/Functor.h"
#include "Math/Integrator.h"
#include "Math/WrappedMultiTF1.h"
#include "Math/WrappedParamFunction.h"
#include "Math/GSLMCIntegrator.h"
#include <ROOT/RDataFrame.hxx>
#include "TStopwatch.h"

auto sqr(double x)
{ return x * x; }

double f_1(double x) {
Double_t par={0.504536,0.0671194,6829.58};
Double_t val = (sqr(par)*par)/sqr(par-par)*(exp(-par*x)*((par-par)*x-1)+exp(-par*x))*par;
return val;
}

double fint_1()
{
ROOT::Math::Functor1D wf(&f_1);
ig.SetFunction(wf);
Double_t T=60.0;
Double_t fintval = ig.Integral(0.0,T);

return fintval;
}

double f(double x) {
Double_t p0=0.504536, p1=0.0671194, p2=6829.58, p3=fint_1();
return ((sqr(p0)*p1/sqr(p1-p0))*(exp(-p0*x)*((p1-p0)*x-1.0)+exp(-p1*x))*p2)/p3;
}
//____________________________

double fint1(double t)
{
ROOT::Math::Functor1D wf(&f);
ig.SetFunction(wf);
Double_t tau=0.9;
Double_t fintval = ig.Integral(t,t+tau);

return fintval;
}
//_______________________________i=6-------------------------------------------------------

double ff6(const double *x, const double *par)
{
const int i=6;
Double_t tau=0.9;
Double_t T=60.0;
Double_t t2min=x+tau;
Double_t t2max=T-(i-2.0)*tau;
Double_t t2=(t2max-t2min)*x+t2min;
Double_t dt2=t2max-t2min;
Double_t t3min=t2+tau;
Double_t t3max=T-(i-3.0)*tau;
Double_t t3=(t3max-t3min)*x+t3min;
Double_t dt3=t3max-t3min;
Double_t t4min=t3+tau;
Double_t t4max=T-(i-4.0)*tau;
Double_t t4=(t4max-t4min)*x+t4min;
Double_t dt4=t4max-t4min;
Double_t t5min=t4+tau;
Double_t t5max=T-(i-5.0)*tau;
Double_t t5=(t5max-t5min)*x+t5min;
Double_t dt5=t5max-t5min;
Double_t t6min=t5+tau;
Double_t t6max=T-(i-6.0)*tau;
Double_t t6=(t6max-t6min)*x+t6min;
Double_t dt6=t6max-t6min;
Double_t ffval = (TMath::Factorial(par)/TMath::Factorial(par-i)) *f(x) *f(t2) *f(t3) *f(t4) *f(t5) *f(t6) *dt2 *dt3 *dt4 *dt5 *dt6 *pow(fint1(x) +fint1(t2) +fint1(t3) +fint1(t4) +fint1(t5) +fint1(t6),par-i);
return ffval;
}
double S6j(double j)
{
const int i=6;
Double_t tau=0.9;
Double_t T=60.0;
Double_t epsabs=1e-2;
ROOT::Math::WrappedParamFunction<> wf( &ff6, i, 1);
Double_t p = {j};
wf.SetParameters(p);
// Create the Integrator
//ROOT::Math::GSLMCIntegrator ig(ROOT::Math::MCIntegration::kVEGAS);
// Set parameters of the integration
ig.SetFunction(wf);
ig.SetAbsTolerance(epsabs);
Double_t xmin[i] = {0.0,0.0,0.0,0.0,0.0,0.0};
Double_t xmax[i] = {T-(i-1.0)*tau,1.0,1.0,1.0,1.0,1.0};
Double_t Sjvalue =ig.Integral(xmin, xmax);
cout << "Error: " << ig.Error() << endl;
cout << "Matrix element S66: ";
return Sjvalue;
}
TStopwatch timer;
void spencerIntegrals_rootforumtest() {
ROOT::EnableImplicitMT(0);
timer.Start();
cout << S6j(6) << endl;
timer.Stop();
printf("---RT=%7.5f m, RT=%7.5f h\n",timer.RealTime()/60.0,timer.RealTime()/60.0/60.0);
}
``````

Hello Jeremias,

First of all your integration code contains a large overhead by creating the `ROOT::Math::Integrator` class inside the inner function. Since you are always integrating the same function, the class should be created outside the inner computations and this will speed up by order of magnitude your integral.
I attach the speed-up code. Vegas takes longer probably because it requires more function call in your case. You can probably tune a bit by using different parameters, like a reduced number of iterations, since you require a rather large tolerance. In the attached code there is an example how to tune the parameters.
If you have further questions, please let me know

Best regards

Lorenzo

testMC.C (3.9 KB)

1 Like

Thank you @moneta! It works now, it is much faster, and after increasing the maximum number of calls to around 1E6 I was able to get a two orders of magnitude better uncertainty in the result.

Cheers,
Jeremias

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.