Benchmarking expressions

New method

Mathematica 10 introduced a Benchmark and BenchmarkPlot functions in the included package GeneralUtilities. The latter automates and extends the process described in the Legacy section below. Version 10.1 also introduced RepeatedTiming which is like a refined version of timeAvg that uses TrimmedMean for more accurate results. Here is an overview of these functions.

Before using BenchmarkPlot one must load the package with Needs or Get. Its usage message reads:

Benchmark has a similar syntax but produces a table of numeric results rather than a plot.

The Options and default values for BenchmarkPlot are:

{
 TimeConstraint -> 1.`,
 MaxIterations -> 1024, 
 "IncludeFits" -> False,
 "Models" -> Automatic
}

Notes:

• BenchmarkPlot uses caching therefore if a function is changed after benchmarking its timings may not be accurate in future benchmarking. To clear the cache one may use:

  Clear[GeneralUtilities`Benchmarking`PackagePrivate`$TimingCaches]
  

• There is a bug in 10.1.0 regarding BenchmarkPlot; a workaround is available.

Legacy method

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.

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 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.