Walking Randomly

Mike, so you are saying that fitting 7 datapoints takes 0.5 seconds ?! That sounds ridiculous right? I tried doing it in Mathematica, it takes 0.1 second the first time, but after the packages are loaded only milliseconds. And that even includes 2 plots!

AbsoluteTiming[
x={0,0.2,0.4,0.6,0.75,0.9,1};
y={1,1.22140275816017,1.49182469764127,1.822118800390509,2.117000016612675,2.45960311115695,2.718281828459045};
l1=ListPlot[Transpose[{x,y}]];
timesol=Interpolation[Transpose[{x,y}],InterpolationOrder->1, Method->”Hermite”];
p1=Plot[timesol[x],{x,0,1}];
Show[{l1,p1}]
]

So maybe she should give Mathematica a try? What do you think?

Hi Mike,

I enjoy reading your blog.

Your timing of INTERP1 of Matlab surprises me. I assume you use the latest version of Matlab.
I use R2008a on a 6+ years old Dell Precission 450 with Windows XP SP3.

I think that “… individual timings can vary wildy …” is not true regarding Matlab. (Assuming that not other processes are competing for the cpu-cycles.)

Here are a couple of my results.

Execution of your original script:

>> clear all
…
Elapsed time is 0.260425 seconds. — first call – functions are loaded into memory
Elapsed time is 0.048815 seconds. — second call – ???
Elapsed time is 0.011384 seconds.
Elapsed time is 0.009701 seconds.
Elapsed time is 0.009132 seconds.
Elapsed time is 0.009301 seconds.
Elapsed time is 0.008757 seconds.

I repeated this a few time and see the same pattern every time.

Execution two times of the loop. It does not include the assignment of t. Mean and standard deviation of the results of the second run are:

>> mean(matlabtimes)
ans =
0.0082
>> std(matlabtimes)
ans =
6.7250e-004

/ per isakson

Mike,

I didn’t see that, sorry for the quick judging. However don’t use table in Mathematica, it is most often not so fast, there are faster ways:
timesol /@ Range[0, 1, 0.00005]
or even
timesol[Range[0,1,0.00005]] (as the function automatically threads over a list of arguments).

Please try those as well, it could be quite a bit of difference!, and while your at it; please also average over 100 runs ;-)

Sander

Mike,

I often try to circumvent using table, but rather using Map (/@). Note that Table is not always slower, it is sometimes faster than Map. But if you have a list of values, and you want to evaluate a function at those values, most often Map (/@) is the fastest. It just depends on the task which is faster, I often try both with varying result.

I think in general you can say: For < Table and in some cases For < Table < Map.

This code finds the first n primes as a list, it can be done in many, many ways, here are a couple:

n=600000;
t1=AbsoluteTiming[Prime[Range[1,n]];];
t2=AbsoluteTiming[Prime/@Range[1,n];];
t3=AbsoluteTiming[Map[Prime,Range[1,n]];];
t4=AbsoluteTiming[Table[Prime[x],{x,1,n}];];
t5=AbsoluteTiming[Reap[Do[Sow[Prime[i]],{i,1,n}]][[2,1]];];
t6=AbsoluteTiming[res=ConstantArray[0,n];Do[res[[i]]=Prime[i],{i,1,n}];];
t7=AbsoluteTiming[Reap[For[i=1,i<=n,i++,Sow[Prime[i]]]][[2,1]];];
t8=AbsoluteTiming[res=ConstantArray[0,n];For[i=1,i<=n,i++,res[[i]]=Prime[i]];];
times={t1,t2,t3,t4,t5,t6,t7,t8}[[All,1]]
BarChart[times]

Sometimes it is worth looking in to the way you make a list!

Sander

Counterexample when Table is faster than Map could be the simple task of adding 2:

AbsoluteTiming[# + 2 & /@ Range[10^6];]
AbsoluteTiming[Table[i + 2, {i, 10^6}];]

(note that Range[10^6)+2 is way faster, or one could even type Range[3,10^6+2], which is again much faster)

Cheers,

Sander

An example where Map is faster than Table would be

AbsoluteTiming[Table[PrimeQ[i], {i, 0, 10000000}];]
AbsoluteTiming[PrimeQ /@ Range[0, 10000000];]

On my computer the latter is about 7.5% faster. Not a lot, but it’s something :-)

Ok, last comment of the day!

Sander

Mike,

Now, I have run the loop with the larger t (0.2e6 elements, numel(t)=200001 ). It takes 10 times longer. In this case, function or script doesn’t affect the speed – I checked. (The JIT-compiler used to be more effective with functions.)

>> clear all
mean = 0.0858, std = 0.0257
mean = 0.0824, std = 0.0015

/ per

I never trust numerical “benchmark” demonstrations that take less than a second to run, much less those that take a fraction of a second. There could be too much stuff going on in the background. If you want to benchmark something, write a loop that does the same thing over and over, close all applications, let it run for several minutes.

I do believe that NAG should be much faster. A lot of Matlab functions are written in Matlab’s language, but NAG library is compiled code. Now, if only my university had a license for NAG. It would be interesting to see how much faster is NAG’s optimization code, arma estimation, and other statistical routines. I suspect much fast since most of this code in Matlab’s toolboxes is based on interpreted m-files. This is specially bad for optimization routines, most of which have to run in a loop.

Oh, I am sorry. I did not intend to ask you to benchmark anything. I just wondered, in general, how much, if at all, NAG functions could be faster than the MATLAB implementations for doing statistical work in things like say time series analysis. I am almost certain that NAG will be faster since probably all of NAG is compiled code but only some of MATLAB’s functions are compiled into libraries. For example, Matlab’s ARMA-GARCH family of functions in econometrics toolbox seems to be written entirely as MATLAB scripts, as well as many things in the optimization toolbox. I believe there is a general view in the statistics community that MATLAB’s optimization is weaker than other packages, but I don’t know to what extent this is true. When I have time, I’ll download the trial copy of NAG toolbox to try it myself.