15
$\begingroup$

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?

$\endgroup$
2
  • 2
    $\begingroup$ @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. $\endgroup$ Mar 31, 2013 at 20:47
  • $\begingroup$ @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 :) $\endgroup$
    – rm -rf
    Mar 31, 2013 at 23:00

3 Answers 3

19
$\begingroup$

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:

enter image description here

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]

Mathematica graphics

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

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

ListLinePlot[timings2]

Mathematica graphics

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.

$\endgroup$
6
  • 2
    $\begingroup$ 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.) $\endgroup$
    – Mr.Wizard
    Mar 31, 2013 at 21:33
  • $\begingroup$ thanks for this answer and congratulations on milestone...thanks for all instruction/corrections/editing/moderation and vigilance $\endgroup$
    – ubpdqn
    May 14, 2015 at 1:07
  • $\begingroup$ @Mr.Wizard Why I can't find BenchmarkPlot in 10.2? $\endgroup$
    – matheorem
    Oct 11, 2015 at 12:10
  • $\begingroup$ @matheorem You have loaded the package with Needs["GeneralUtilities`"]? I would be surprised if they removed it, but I do not have version 10.2 myself so I cannot check. $\endgroup$
    – Mr.Wizard
    Oct 12, 2015 at 7:00
  • $\begingroup$ @Mr.Wizard Thank you, about timeAvg, I have a question. What if when running the second time and after, the expression uses cache? Because according to the doc, we'd better use ClearSystemCache[] before every timing. To circumvent this problem, I don't know if clear cache in every loop is good, for clear cache may also consumes considerable time. Then I can only think of this: to prepare enough data set before timing, and run through these data every fresh loop. Looking forward to your opinion. $\endgroup$
    – matheorem
    Oct 13, 2015 at 14:16
4
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ Don't forget AbsoluteTiming. $\endgroup$ Mar 31, 2013 at 20:44
4
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.