random variate from power exponential distribution: continue

Mon Oct 18 16:43:00 GMT 2004

Hi,
again, let me apologize for the length of this email. This email is
about the generation of power exponential pseudo random variates and its
implementation inside GSL. If you are not interested in this issue, you
can safely ignore it:-). For the rest of you, here's a table of
contents:

Sec.1 Good news: the code can be rewritten to be much faster

Sec.2 Bad news: the relative performances of the methods seem to depend
on the platform

Section 1
---------

As I wrote in a previous email, the function for the generation of power
exponential random variables actually implemented inside GSL (file
randist/exppow.c) can be improved by exploiting the analytic expression
of the threshold levels for the rejection method, as provided in
Tadikamalla's paper (see my previous email for details). Indeed,
rewriting the code, I obtained a noticeable increase in speed. Moreover,
when I tested them on my laptop, the three different methods, namely the
transformation of a gamma variate(GS), the rejection method from a
Laplace distribution (ED) and the rejection method from a Gaussian (EN),
performed as described by Tadikamalla (and as implicitly assumed inside
the actual GSL implementation). Here the times for the generation of 1
million variates as a function of b using the new version of the three
methods(ED can only be used for b>=1 and EN for b>=2)

#b      GS      ED      EN
0.500   0.560   nan     nan
0.700   1.600   nan     nan
0.900   1.500   nan     nan
1.001   1.340   1.140   nan
1.250   1.350   1.120   nan
1.500   1.340   1.180   nan
1.750   1.340   1.230   nan
2.001   1.310   1.200   1.100
2.250   1.300   1.240   1.130
2.500   1.290   1.270   1.180
2.750   1.270   1.300   1.200
3.000   1.270   1.150   1.080
3.500   1.260   1.380   1.280
4.000   1.240   1.330   1.260
5.000   1.220   1.390   1.320
6.000   1.200   1.430   1.370
8.000   1.170   1.530   1.430

as can be seen the optimal method as a function of b is:

    b<1       GS
1<= b < 2     ED
2<= b <~ 3.5  EN
    B >~ 3.5  GS

broadly in line with Tadikamalla's original findings.

Section 2
---------

The previous results suggest for a simple rewriting of the present GSL
code. However, I decided to benchmark the routines also on different
machines. What I found surprised me: on AMD systems, the relative
performances of the tree methods seem different. Here the results of my
test program on an Athlon XP

#b      GS      ED      EN
0.500   0.560   nan     nan
0.700   1.370   nan     nan
0.900   1.250   nan     nan
1.001   1.180   1.040   nan
1.250   1.170   1.110   nan
1.500   1.170   1.140   nan
1.750   1.160   1.200   nan
2.001   1.140   1.220   1.170
2.250   1.130   1.270   1.180
2.500   1.120   1.290   1.210
2.750   1.100   1.330   1.240
3.000   1.090   1.190   1.140
3.500   1.080   1.430   1.340
4.000   1.060   1.240   1.190
5.000   1.050   1.310   1.250
6.000   1.040   1.390   1.320
8.000   1.000   1.440   1.370

and on a dual opteron

#b      GS      ED      EN
0.500   0.300   nan     nan
0.700   1.010   nan     nan
0.900   0.910   nan     nan
1.001   0.840   0.770   nan
1.250   0.850   0.870   nan
1.500   0.830   0.910   nan
1.750   0.820   0.950   nan
2.001   0.810   0.970   0.880
2.250   0.850   1.150   1.010
2.500   0.850   1.180   1.050
2.750   0.830   1.150   0.960
3.000   0.770   1.120   0.980
3.500   0.760   1.150   1.020
4.000   0.750   1.200   1.050
5.000   0.730   1.280   1.110
6.000   0.730   1.300   1.170
8.000   0.710   1.410   1.250

In the first case method EN is always suboptimal; in the second case
method GS always runs faster.

Notice that these results have been obtained running EXACTLY the same
source code, that you can find in attachment. The source has always been
compiled with

gcc -O2 test_exppow.c -lgsl -lgslcblas -o test

Moreover, on all systems, the GSL were installed with the standard
configure-make-make install sequence (i.e. CFLAGS = "-g -O2").

I hope that my surprise about this finding is due to my ignorance and
that someone could easily explain me the reason of these differences
across systems. It would be nice if someone can run the attached program
on his/her machines (just compile it with the line above and issue the
command "./test" without options) and let me know about the results. In
any case, I'll prepare a replacement for exppow.c with the
improvements mentioned in Section 1. These should lead to FASTER
programs on any system.

Best,
	Giulio.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: test_exppow.c
Type: text/x-csrc
Size: 3668 bytes
Desc: not available
URL: <http://sourceware.org/pipermail/gsl-discuss/attachments/20041018/c34eeadc/attachment.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: not available
URL: <http://sourceware.org/pipermail/gsl-discuss/attachments/20041018/c34eeadc/attachment.sig>