# Benchmarking expressions

It often happens that I'd like to know which of two or more alternative expressions (that should evaluate to the same value) is fastest. For example, below are three ways to define a function:

fooVersion1[x_]/; x > 0 := E^(-1/x);
fooVersion1[x_] := 0;

fooVersion2[x_] := If[x > 0, E^(-1/x), 0];

fooVersion3 = Function[x, If[x > 0, E^(-1/x), 0]];


(Granted, strictly speaking, these three versions are not completely equivalent, but at least, for all real scalars x, fooVersion1[x], fooVersion2[x], and fooVersion3[x] should evaluate to the same value.)

Is there a way to determine which one is fastest that is both simple and reasonably reliable?

-
@rm-rf Not sure about the dupe. Benchmarking and profiling are certainly very related, yet I would hesitate to call them the same thing. Benchmarking is usually a process of comparing several different solutions in terms of speed, while profiling is more about distribution of running time for a single function over its sub-function calls. –  Leonid Shifrin Mar 31 '13 at 20:47
@LeonidShifrin Ok, fair enough; I agree with you. If it eventually gets closed, please let me know and I'll reopen it. If not, the votes will expire away after a few days :) –  rm -rf Mar 31 '13 at 23:00

I described my basic method of comparative benchmarking in this answer:
How to use "Drop" function to drop matrix' rows and columns in an arbitrary way?

I shall again summarize it here. I make use of my own modification of Timo's timeAvg function which I first saw here. It is:

SetAttributes[timeAvg, HoldFirst]

timeAvg[func_] := Do[If[# > 0.3, Return[#/5^i]] & @@ Timing@Do[func, {5^i}], {i, 0, 15}]


This performs enough repetitions of a given operation to exceed the threshold (0.3 seconds) with the aim getting a sufficiently stable and precise timing.

I then use this function on input of several different forms, to try to get a first order mapping of the behavior of the function over its domain. If the function accepts vectors or tensors I will try it with both packed arrays and unpacked (often unpackable) data as there can be great differences between methods depending on the internal optimizations that are triggered.

I build a matrix of timings for each data type and plot the results:

funcs = {fooVersion1, fooVersion2, fooVersion3};

dat = 20^Range[0, 30];

timings = Table[timeAvg[f @ d], {f, funcs}, {d, dat}];

ListLinePlot[timings]


Here with Real data as a second type (note the N@):

timings2 = Table[timeAvg[f @ d], {f, funcs}, {d, N@dat}];

ListLinePlot[timings2]


One can see there is little relation between the size of the input number and the speed of the operation with either Integer or machine-precision Real data.

The sawtooth shape of the plot lines is likely the result of insufficient precision (quantization). It could be reduced by increasing the threshold in the timeAvg function, or better in this case to time multiple applications of each function under test in a single timeAvg operation.

-
I just noticed this is my one thousandth answer on Mathematica.SE. :-) (Some users will see less listed in my profile as it does not show them deleted answers.) –  Mr.Wizard Mar 31 '13 at 21:33

The Timing function does this kind of thing nicely. For instance:

Timing[Table[fooVersion1[RandomReal[{-1, 1}]], {i, 1, 1000000}];]
Timing[Table[fooVersion2[RandomReal[{-1, 1}]], {i, 1, 1000000}];]
Timing[Table[fooVersion3[RandomReal[{-1, 1}]], {i, 1, 1000000}];]


which shows that the first is a bit faster than the second which is a bit faster than the third. I got 2.8, 3.4, and 4.2 seconds.

-
Don't forget AbsoluteTiming. –  Sjoerd C. de Vries Mar 31 '13 at 20:44

Don't forget the Which and Piecewise versions:

fooVersion4[x_] = Which[x > 0, E^(-1/x), True, 0]

fooVersion5[x_] = Piecewise[{{E^(-1/x), x > 0}}]


Timings on my computer:

AbsoluteTiming[Table[fooVersion1[RandomReal[{-1, 1}]], {i, 1, 1000000}];]
AbsoluteTiming[Table[fooVersion2[RandomReal[{-1, 1}]], {i, 1, 1000000}];]
AbsoluteTiming[Table[fooVersion3[RandomReal[{-1, 1}]], {i, 1, 1000000}];]
AbsoluteTiming[Table[fooVersion4[RandomReal[{-1, 1}]], {i, 1, 1000000}];]
AbsoluteTiming[Table[fooVersion5[RandomReal[{-1, 1}]], {i, 1, 1000000}];]
{3.330165, Null}
{4.132347, Null}
{4.970121, Null}
{4.008108, Null}
{5.052208, Null}


Unsurprisingly, Piecewise is the slowest. To my surprise, Which is slightly faster than If.

I would probably try to avoid the first type of input ... fooVersion1 ... because Mathematica may have trouble operating on it symbolically. For example, compare the correct result:

Integrate[fooVersion2[x], {x, -Infinity, 1}]
1/E + ExpIntegralEi[-1]


... to what happens when you use fooVersion1:

Integrate[fooVersion1[x], {x, -Infinity, 1}]
0


... and saving a second once a year (and getting the wrong result) is not worth the risk.

-