while profiling a simulation application that makes use of itpp, i've noticed that the itpp::Normal_RNG::sample function is showing up fairly high in the profile rankings (~5%), along with some of the child functions it calls, like log() from libm. while looking into this i came across an article about an alternative normal RNG algorithm called the "ziggurat method", devolped by marsaglia and later tsang, which is supposed to be noticably faster than the variants of box-mueller. indeed, it seems that at some point in the distant past the folks at matlab realized this and now they too are using this algorithm for their normal RNG. i don't have a copy of TAoCP handy, but i've heard there's some reference to this there as well.
Sure! If it is confirmed to work for Matlab and the speed improvement will be noticeable, we will for sure consider replacing our current implementation with this "ziggurat method"...
BR,
/ediap
PS. I am quite busy this week, so I will probably go back to this issue in the nearest future. Thanks for your patience.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
which is GPL'd, even (some of the other code was non-commercial only).
i've incorporated this into the Normal_RNG, replacing the Box-Mueller algorithm. the patch is fairly clean, changing the guts of sample and a few new private const/static tables (holding various precomputed variables related to the ziggurat steps, i defer to the algorithm here).
running a simple benchmarking application that sums 100,000,000 samples from the Normal_RNG takes just a bit shy of 12 seconds
with the original algorithm, and with the ziggurat method it is around 2-3 times faster, averaging around 4.5 seconds.
these tests were done with
- the latest development release (i don't have admin rights and couldn't get the various autotools stuff setup to build from svn)
- g++ 4.1.2 with -O3
- RHEL 3.x running a 2.4 kernel
- a dedicated 3Ghz Xeon processor
though i would expect to see such speedups with other OS/compiler combinations as well. where should i send a patch for review?
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
okay, i'll pretty up the patch with some comments and post it to the tracker.
before going any further, i should give the disclaimer that i'm not a mathematician/statistician, so you should probably have somebody with a better knowledge of such thigns audit the ideas and resulting code that i present. also, do you have any tests that you use for your RNG's to ensure the distribution qualities?
now, back to the matter at hand wrt complex normal distribution, according to wikipedia[1] you can generate an n-variate normal distributed vector by:
Compute the Cholesky decomposition (matrix square root) of Sigma (the covariance matrix), that is, find the unique lower triangular matrix A such that AA^T = Sigma
Let Z=(z1,...,zn)^T be a vector whose components are n independent standard normal variates (which can be generated, for example, by using the Box-Muller transform).
Let X be mu (the mean) + AZ.
in our case, we have two dimensions and a single variance, so Sigma would be a 2x2 diagonal matrix with diagonals equal to the variance.
A would then be a diagonal matrix of diagonals sigma = std::sqrt(variance).
Z could simply be two variants generated from a Normal_RNG.
mu is represented by as a complex number via the constructor.
then we have:
a = normal_rng.sample()
b = normal_rng.sample()
X = std::complex( (mu.real() + sigmaa), (mu.imag() + sigmab) )
In order to generate a complex normal distribution with independent real and imaginary parts (diagonal covariance matrix SIGMA), it is really as simple as seanius wrote:
a = normal_rng.sample()
b = normal_rng.sample()
X = std::complex( (mu.real() + sigmaa), (mu.imag() + sigmab) )
In this case, both real and imaginary parts have the same variance sigma^2.
In my humble opinion, the itpp function should generate a zero mean (mu=0) and unitary variance (sigma=1) complex normal, as in
a = normal_rng.sample()
b = normal_rng.sample()
X = std::complex( a, b )
But, once more, it's just my 2 cents...
BR,
daniloz
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
first thanks for the comments daniloz. wrt the mean/variance, i think that the sample() functions return just this, and that then the public accessor functions like operator() return this number modified according to the mean/variance provided to the RNG constructor (this was the cause my confusion in the previous postings).
second, to adam: you want a patch for this too?
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
It would be nice to verify the performance difference between the current implementation of the Complex_Normal_RNG and the new one using the Ziggurat based Normal_RNG. And of course, patches are always welcome, since they simplify our (my ;-) maintenance work.
BR,
/ediap
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
i went ahead and did this, but my implementation has raised a quick question that i'd like a second opinion on just to make sure i'm understanding things correctly (again insert previous disclaimer).
for the variance of a complex normal distribution where mean=0, variance=1, should variance=1 for the series of real and imaginary numbers independantly, or should it be for the complex number series as a whole?
with the existing generator, the individual members have mean0, variance=~0.5, but if i calculate the variance as
E[ (X-mean)(X-mean)* ] where foo* is the complex conjugate
i get a variance of 1, which is what i guess should be expected.
using the simple approach mentioned earlier in this thread:
return complex(normal.sample(),normal.sample())
i get mean=0, variance = 2 for the whole series, but variance = 1 for the individual members.
however, if i divide both parts by sqrt(2), the variance seems to match the existing generator. so then, i guess the following would
suffice for a complex normal series?
where normal is a Normal_RNG which replaces the standard private RNG member in the class, and oneoversq2 is 1/sqrt(2), which could be calculated once in the constructor (or just hard-coded).
does this seem correct?
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
You are right with your assumptions. The variance should describe the complex number, not its real and imaginary parts separately. So you need to scale each Normal_RNG part by 1.0/sqrt(2.0) factor.
BR,
/ediap
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
okay, i've created a new ticket in the tracker for this, attaching a patch and also the benchmarking program i've used with some sample timings.
another question (i seem to keep finding more): all of these rng's are seperate and unrelated classes... any reason for the lack of oo-inheritance and/or templating? the latter would of course end up in a non-backwards-compatible api change, but the former could be used at the least to
reduce the LoC a bit.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
The answer to your question, especially the suggestion to use some kind of inheritance, is simple: performance.
Some time ago I tried to redesign the Random_Generator class, so it could have a common interface (base abstract class) with a set of implementations in its descendent classess (classical Mersenne-Twister, SFMT, etc.). However, the performance loss due to virtual functions and pointers were quite high (almost four time slower generation that the present code). So, I decided to keep this code as much inlined and independent as possible.
BR,
/ediap
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
hey again,
while profiling a simulation application that makes use of itpp, i've noticed that the itpp::Normal_RNG::sample function is showing up fairly high in the profile rankings (~5%), along with some of the child functions it calls, like log() from libm. while looking into this i came across an article about an alternative normal RNG algorithm called the "ziggurat method", devolped by marsaglia and later tsang, which is supposed to be noticably faster than the variants of box-mueller. indeed, it seems that at some point in the distant past the folks at matlab realized this and now they too are using this algorithm for their normal RNG. i don't have a copy of TAoCP handy, but i've heard there's some reference to this there as well.
some links:
http://www.jstatsoft.org/v05/i08/ziggurat.pdf
http://www.mathworks.com/company/newsletters/news_notes/clevescorner/spring01_cleve.html
http://en.wikipedia.org/wiki/Ziggurat_algorithm
would you be interested in such an alternative normal RNG?
Sure! If it is confirmed to work for Matlab and the speed improvement will be noticeable, we will for sure consider replacing our current implementation with this "ziggurat method"...
BR,
/ediap
PS. I am quite busy this week, so I will probably go back to this issue in the nearest future. Thanks for your patience.
after spending the morning reading a few papers and digging around for sample code, i've managed to implement this.
i found some example c code here:
http://seehuhn.de/comp/ziggurat/
which is GPL'd, even (some of the other code was non-commercial only).
i've incorporated this into the Normal_RNG, replacing the Box-Mueller algorithm. the patch is fairly clean, changing the guts of sample and a few new private const/static tables (holding various precomputed variables related to the ziggurat steps, i defer to the algorithm here).
running a simple benchmarking application that sums 100,000,000 samples from the Normal_RNG takes just a bit shy of 12 seconds
with the original algorithm, and with the ziggurat method it is around 2-3 times faster, averaging around 4.5 seconds.
these tests were done with
- the latest development release (i don't have admin rights and couldn't get the various autotools stuff setup to build from svn)
- g++ 4.1.2 with -O3
- RHEL 3.x running a 2.4 kernel
- a dedicated 3Ghz Xeon processor
though i would expect to see such speedups with other OS/compiler combinations as well. where should i send a patch for review?
Please open a new request at Feature Request Tracker and attach the patch with a short description there.
BTW, have you found something related to the generation of a complex normal distribution? Maybe some better algorithm exists as well?
BR,
/ediap
okay, i'll pretty up the patch with some comments and post it to the tracker.
before going any further, i should give the disclaimer that i'm not a mathematician/statistician, so you should probably have somebody with a better knowledge of such thigns audit the ideas and resulting code that i present. also, do you have any tests that you use for your RNG's to ensure the distribution qualities?
now, back to the matter at hand wrt complex normal distribution, according to wikipedia[1] you can generate an n-variate normal distributed vector by:
in our case, we have two dimensions and a single variance, so Sigma would be a 2x2 diagonal matrix with diagonals equal to the variance.
A would then be a diagonal matrix of diagonals sigma = std::sqrt(variance).
Z could simply be two variants generated from a Normal_RNG.
mu is represented by as a complex number via the constructor.
then we have:
a = normal_rng.sample()
b = normal_rng.sample()
X = std::complex( (mu.real() + sigmaa), (mu.imag() + sigmab) )
it can't be that simple, can it?
[1] http://en.wikipedia.org/wiki/Multivariate_normal_distribution
oh, i think i've managed to confuse myself between what sample() does and what operator() does in the context of the Complex_Normal_RNG.
so if sample() gives the raw values, couldn't we just return a complex number generated from two Normal_RNG numbers?
Hello all,
In order to generate a complex normal distribution with independent real and imaginary parts (diagonal covariance matrix SIGMA), it is really as simple as seanius wrote:
a = normal_rng.sample()
b = normal_rng.sample()
X = std::complex( (mu.real() + sigmaa), (mu.imag() + sigmab) )
In this case, both real and imaginary parts have the same variance sigma^2.
In my humble opinion, the itpp function should generate a zero mean (mu=0) and unitary variance (sigma=1) complex normal, as in
a = normal_rng.sample()
b = normal_rng.sample()
X = std::complex( a, b )
But, once more, it's just my 2 cents...
BR,
daniloz
hi guys,
first thanks for the comments daniloz. wrt the mean/variance, i think that the sample() functions return just this, and that then the public accessor functions like operator() return this number modified according to the mean/variance provided to the RNG constructor (this was the cause my confusion in the previous postings).
second, to adam: you want a patch for this too?
> second, to adam: you want a patch for this too?
It would be nice to verify the performance difference between the current implementation of the Complex_Normal_RNG and the new one using the Ziggurat based Normal_RNG. And of course, patches are always welcome, since they simplify our (my ;-) maintenance work.
BR,
/ediap
hey adam,
i went ahead and did this, but my implementation has raised a quick question that i'd like a second opinion on just to make sure i'm understanding things correctly (again insert previous disclaimer).
for the variance of a complex normal distribution where mean=0, variance=1, should variance=1 for the series of real and imaginary numbers independantly, or should it be for the complex number series as a whole?
with the existing generator, the individual members have mean0, variance=~0.5, but if i calculate the variance as
E[ (X-mean)(X-mean)* ] where foo* is the complex conjugate
i get a variance of 1, which is what i guess should be expected.
using the simple approach mentioned earlier in this thread:
return complex(normal.sample(),normal.sample())
i get mean=0, variance = 2 for the whole series, but variance = 1 for the individual members.
however, if i divide both parts by sqrt(2), the variance seems to match the existing generator. so then, i guess the following would
suffice for a complex normal series?
return complex( normal.sample() * oneoversq2 , normal.sample() * oneoversq2 )
where normal is a Normal_RNG which replaces the standard private RNG member in the class, and oneoversq2 is 1/sqrt(2), which could be calculated once in the constructor (or just hard-coded).
does this seem correct?
Hi Sean,
You are right with your assumptions. The variance should describe the complex number, not its real and imaginary parts separately. So you need to scale each Normal_RNG part by 1.0/sqrt(2.0) factor.
BR,
/ediap
okay, i've created a new ticket in the tracker for this, attaching a patch and also the benchmarking program i've used with some sample timings.
another question (i seem to keep finding more): all of these rng's are seperate and unrelated classes... any reason for the lack of oo-inheritance and/or templating? the latter would of course end up in a non-backwards-compatible api change, but the former could be used at the least to
reduce the LoC a bit.
The answer to your question, especially the suggestion to use some kind of inheritance, is simple: performance.
Some time ago I tried to redesign the Random_Generator class, so it could have a common interface (base abstract class) with a set of implementations in its descendent classess (classical Mersenne-Twister, SFMT, etc.). However, the performance loss due to virtual functions and pointers were quite high (almost four time slower generation that the present code). So, I decided to keep this code as much inlined and independent as possible.
BR,
/ediap