[Gretl-devel] simulation speed
Hi,
yet another useless benchmarking exercise: I was playing around a little
comparing the speed of various matrix languages with respect to doing
stochastic simulations. The codes below draw repeated random samples,
compute their means, and find the 95% percent quantile of the empirical
distribution of the means. I distinguish between a "naive" loop
programming variant and one working only with matrix functions. Here are
the results on Ubuntu Linux 9.10 (I ran the scripts/programs several
times and they are pretty robust):
Octave 3.0.5:
loop ca. 3.3 seconds
vectorized ca. 0.8 seconds
Python 2.6.4/Numpy 1.3.0/matplotlib 0.99(?):
loop ca. 2.1 seconds
vectorized ca. 1.5 seconds
gretl 1.8.6cvs:
loop ca. 2.6 seconds
vectorized ca. 2.0 seconds
Ox 5.1:
loop ca. 1.1 seconds
vectorized ca. 0.7 seconds
Ox confirms its reputation about speed, but for me it was surprising
that Octave is almost as fast with the vectorized code. (But Octave has
the slowest loops.) Gretl seems to have the smallest relative gain going
from loops to vectorized code. But considering that gretl is quite new
in the field of matrix programming I think it is doing fine.
cheers,
sven
<Octave-code>
length = 1000;
iterations = 10000;
whichp = 0.95;
tic;
# looping variant
averages = [];
for i=1:iterations;
mrandom = randn(length,1);
averages = [averages ; mean(mrandom)];
endfor;
myquant = quantile(averages,whichp);
toc;
printf ("Result: %f\n", myquant);
tic;
# vectorized variant
mrandom2 = randn(length,iterations);
averages2 = mean(mrandom2); # row vector
myquant2 = quantile(averages2',whichp);
toc;
printf ("Result: %f\n", myquant2);
</Octave-code>
<Python-code>
from numpy import matlib as nm
from matplotlib import mlab as ml
import time
length = 1000
iterations = 10000
whichp = 0.95
starttime = time.time()
# looping variant
averages = nm.empty((0,1))
for i in range(iterations):
mrandom = nm.randn(length,1)
averages = nm.vstack((averages, nm.mean(mrandom)))
myquant = ml.prctile(averages,whichp*100) # per cents
endtime1 = time.time()
print "Result: " + str(myquant)
print "Execution time loop variant, length " + str(length) + ",
iterations " + str(iterations) + ": " + str(endtime1-starttime)
# vectorized variant
starttime2 = time.time()
mrandom2 = nm.randn(length,iterations)
averages2 = nm.mean(mrandom2, axis=0)
myquant2 = ml.prctile(averages2,whichp*100)
endtime2 = time.time()
print "Result: " + str(myquant2)
print "Execution time vectorized variant, length " + str(length) + ",
iterations " + str(iterations) + ": " + str(endtime2-starttime2)
</Python-code>
<gretl-code>
length = 1000
iterations = 10000
whichp = 0.95
set echo off
set messages off
set stopwatch
# looping variant (non-vectorized)
matrix averages = {} # will be col vector
loop iterations
matrix mrandom = mnormal(length,1)
matrix averages = averages | meanc(mrandom)
end loop
matrix myquant = quantile(averages, whichp)
time = $stopwatch
printf "Result: %f\n", myquant
printf "Execution time loop variant, length %d, iterations %d: %f\n",
length, iterations, time
time = $stopwatch # don't count printing overhead
# vectorized variant
matrix mrandom2 = mnormal(length,iterations)
matrix averages2 = transp(meanc(mrandom2))
matrix myquant2 = quantile(averages2, whichp)
time = $stopwatch
printf "Result: %f\n", myquant2
printf "Execution time vectorized variant, length %d, iterations %d:
%f\n", length, iterations, time
</gretl-code>
</Ox-code>
#include <oxstd.h>
const decl length = 1000;
const decl iterations = 10000;
const decl whichp = 0.95;
main(){
decl averages, averages2;
decl myquant, myquant2;
decl mrandom, mrandom2;
decl starttime, endtime;
decl i;
starttime = timer();
// looping variant
averages = <>;
for (i=0; i<iterations; ++i)
{ mrandom = rann(length,1);
averages = averages | meanc(mrandom);
}
myquant = quantilec(averages,whichp);
endtime = timer();
print("Result: ", myquant, "\n");
print("elapsed time: ", timespan(starttime,endtime));
starttime = timer();
// vectorized variant
mrandom2 = rann(length,iterations);
averages2 = meanc(mrandom2); // row vector
myquant2 = quantiler(averages2,whichp);
endtime = timer();
print("\nResult: ", myquant2,"\n");
print("elapsed time: ", timespan(starttime,endtime),"\n");
}
</Ox-code>