I am trying to generate random chi-squared numbers in C++, according to some degree of freedom (which can be a float). Several libraries can be used for that purpose, among which the C++11 Standard Library with the header <random>
, or two Boost libraries, boost::math
or boost::random
.
These libraries don't exactly have the same user interface, but more surprisingly, they definitely don't produce numbers at the same speed. I have performed a test with the next toy code:
#include <chrono>
#include <random>
#include <boost/math/distributions/chi_squared.hpp>
#include <boost/random/chi_squared_distribution.hpp>
unsigned int N = 1e6;
double df = 2.4; //degree of freedom for chi-square
//generator and distributions
std::mt19937 gen; //random number engine, boost::mt19937 doesn't change anything
std::chi_squared_distribution<> chi2STL(df); //STL implementation
boost::random::chi_squared_distribution<> chi2BoostRandom(df); //boost::random
boost::math::chi_squared chi2BoostMath(df); //boost::math
boost::uniform_real<> unif(0., 1.); //uniform distribution for boost::math
//taking STL implementation does not change anything
//time measure
auto point1 = std::chrono::high_resolution_clock::now();
//generate chi2 numbers with STL
for(unsigned int i = 0; i < N; ++i)
{
chi2STL(gen);
}
auto point2 = std::chrono::high_resolution_clock::now();
//generate chi2 numbers with boost::random
for(unsigned int i = 0; i < N; ++i)
{
chi2BoostRandom(gen);
}
auto point3 = std::chrono::high_resolution_clock::now();
//generate chi2 numbers with boost::math
for(unsigned int i = 0; i < N; ++i)
{
boost::math::quantile(chi2BoostMath, unif(gen));
}
auto point4 = std::chrono::high_resolution_clock::now();
//get time
std::cerr << "Time STL loop: " << std::chrono::duration_cast<std::chrono::microseconds>(point2 - point1).count()/1000000. << " s" << std::endl;
std::cerr << "Time Boost Random loop: " << std::chrono::duration_cast<std::chrono::microseconds>(point3 - point2).count()/1000000. << " s" << std::endl;
std::cerr << "Time Boost Math loop: " << std::chrono::duration_cast<std::chrono::microseconds>(point4 - point3).count()/1000000. << " s" << std::endl;
With the next output:
Time STL loop: 0.474826 s
Time Boost Random loop: 0.374868 s
Time Boost Math loop: 5.02392 s
And it turns out that the boost::math method is way slower than the two other ones. I investigated a little more this surprising point, using a normal law instead of a chi-square law. In that case all three methods produce numbers at the same speed, which makes me think that the boost::quantile method is not the time-consuming one.
By the way, I checked the produced numbers, which all seem to be consistent with the chi-square law, at least for my needs.
It seems that the method implemented in boost::math is just less efficient than the two other ones. (Sadly, it is the only one for which I found implementation details online: it makes use of numerical calculations of the incomplete gamma function).
Is such a simplistic interpretation correct ? And in that case, shouldn't the boost::math implementation be deprecated ?