24

So Mathematica is different from other dialects of lisp because it blurs the lines between functions and macros. In Mathematica if a user wanted to write a mathematical function they would likely use pattern matching like f[x_]:= x*x instead of f=Function[{x},x*x] though both would return the same result when called with f[x]. My understanding is that the first approach is something equivalent to a lisp macro and in my experience is favored because of the more concise syntax.

So I have two questions, is there a performance difference between executing functions versus the pattern matching/macro approach? Though part of me wouldn't be surprised if functions were actually transformed into some version of macros to allow features like Listable to be implemented.

The reason I care about this question is because of the recent set of questions (1) (2) about trying to catch Mathematica errors in large programs. If most of the computations were defined in terms of Functions, it seems to me that keeping track of the order of evaluation and where the error originated would be easier than trying to catch the error after the input has been rewritten by the successive application of macros/patterns.

6
  • 4
    Mathematica is not a Lisp dialect. Lisp macros also were designed so that they can be fully expanded at compile time and thus have zero runtime overhead. Commented Nov 15, 2010 at 19:29
  • 1
    I think of Mathematica as an M-expression version of lisp, at least its intellectual heritage. I am not sure what to make of your comment though because it was my understanding that macro-expansion was one of the steps in every read-eval-print cycle in lisp implementations. Because of this I am not sure I understand what you mean by no run-time overhead. I included the lisp tag because I was curious if was an issue that had come up in lisp programming discussions in the past. Commented Nov 15, 2010 at 20:48
  • "premature optimization is the root of all evil" :P Commented Nov 15, 2010 at 22:17
  • 3
    Lisp is based on evaluations of expressions. Mathematica is based on term rewriting. Very different approaches. The Lisp macro-expansion is done all at once by a Lisp compiler. Before evaluation. At runtime there are no Lisp expressions to be simplified or rewritten anymore. Commented Nov 16, 2010 at 1:36
  • @RainerJoswig: You could say that Mathematica applies macros and never gets to evaluation in the Lisp sense. Commented Nov 26, 2011 at 16:49

6 Answers 6

19

The way I understand Mathematica is that it is one giant search replace engine. All functions, variables, and other assignments are essentially stored as rules and during evaluation Mathematica goes through this global rule base and applies them until the resulting expression stops changing.

It follows that the fewer times you have to go through the list of rules the faster the evaluation. Looking at what happens using Trace (using gdelfino's function g and h)

In[1]:= Trace@(#*#)&@x Out[1]= {x x,x^2} In[2]:= Trace@g@x Out[2]= {g[x],x x,x^2} In[3]:= Trace@h@x Out[3]= {{h,Function[{x},x x]},Function[{x},x x][x],x x,x^2} 

it becomes clear why anonymous functions are fastest and why using Function introduces additional overhead over a simple SetDelayed. I recommend looking at the introduction of Leonid Shifrin's excellent book, where these concepts are explained in some detail.

I have on occasion constructed a Dispatch table of all the functions I need and manually applied it to my starting expression. This provides a significant speed increase over normal evaluation as none of Mathematica's inbuilt functions need to be matched against my expression.

Sign up to request clarification or add additional context in comments.

9 Comments

+1 ... and thanks for the book too. BTW, help me build up a bibliography on Mathematica answering here programmers.stackexchange.com/questions/19146/…
Dang, I was going to ask that exact question here on SO. I have seen similar questions about several other languages, e.g., stackoverflow.com/q/336715/181759, and unless there is some general rule disfavoring them I think one about Mathematica would be great.
@Timo I've seen discussions about the on-topicity of my q here on SO, now that the Community Wiki tag is gone. I preferred to starve the trolls by posting there. Of course the q is intended to be a reference for all the Mma community here, so feel free to post whatever material interest you. I'll do the housekeeping to maintain a list in the q text to ease searching.
@belisarius You know, I was just about to ask the "Am I just dense, or how do make a question a wiki?" question. Yes, if the wiki tag is gone then a book/manual question makes less sense here, maybe I should move the packages/tools question to programmers as well.
@Timo, your inference that mma is a giant search and replace engine is somewhat justified as when you do something like f[x_] := x_ you get the error RuleDelayed::rhs : Pattern x_ appears on right-hand side of rule f[x_]:>x_. So, at least to some extent, f[x_] is interpreted as a replacement rule.
|
15

My understanding is that the first approach is something equivalent to a lisp macro and in my experience is favored because of the more concise syntax.

Not really. Mathematica is a term rewriter, as are Lisp macros.

So I have two questions, is there a performance difference between executing functions versus the pattern matching/macro approach?

Yes. Note that you are never really "executing functions" in Mathematica. You are just applying rewrite rules to change one expression into another.

Consider mapping the Sqrt function over a packed array of floating point numbers. The fastest solution in Mathematica is to apply the Sqrt function directly to the packed array because it happens to implement exactly what we want and is optimized for this special case:

In[1] := N@Range[100000]; In[2] := Sqrt[xs]; // AbsoluteTiming Out[2] = {0.0060000, Null} 

We might define a global rewrite rule that has terms of the form sqrt[x] rewritten to Sqrt[x] such that the square root will be calculated:

In[3] := Clear[sqrt]; sqrt[x_] := Sqrt[x]; Map[sqrt, xs]; // AbsoluteTiming Out[3] = {0.4800007, Null} 

Note that this is ~100× slower than the previous solution.

Alternatively, we might define a global rewrite rule that replaces the symbol sqrt with a lambda function that invokes Sqrt:

In[4] := Clear[sqrt]; sqrt = Function[{x}, Sqrt[x]]; Map[sqrt, xs]; // AbsoluteTiming Out[4] = {0.0500000, Null} 

Note that this is ~10× faster than the previous solution.

Why? Because the slow second solution is looking up the rewrite rule sqrt[x_] :> Sqrt[x] in the inner loop (for each element of the array) whereas the fast third solution looks up the value Function[...] of the symbol sqrt once and then applies that lambda function repeatedly. In contrast, the fastest first solution is a loop calling sqrt written in C. So searching the global rewrite rules is extremely expensive and term rewriting is expensive.

If so, why is Sqrt ever fast? You might expect a 2× slowdown instead of 10× because we've replaced one lookup for Sqrt with two lookups for sqrt and Sqrt in the inner loop but this is not so because Sqrt has the special status of being a built-in function that will be matched in the core of the Mathematica term rewriter itself rather than via the general-purpose global rewrite table.

Other people have described much smaller performance differences between similar functions. I believe the performance differences in those cases are just minor differences in the exact implementation of Mathematica's internals. The biggest issue with Mathematica is the global rewrite table. In particular, this is where Mathematica diverges from traditional term-level interpreters.

You can learn a lot about Mathematica's performance by writing mini Mathematica implementations. In this case, the above solutions might be compiled to (for example) F#. The array may be created like this:

> let xs = [|1.0..100000.0|];; ... 

The built-in sqrt function can be converted into a closure and given to the map function like this:

> Array.map sqrt xs;; Real: 00:00:00.006, CPU: 00:00:00.015, GC gen0: 0, gen1: 0, gen2: 0 ... 

This takes 6ms just like Sqrt[xs] in Mathematica. But that is to be expected because this code has been JIT compiled down to machine code by .NET for fast evaluation.

Looking up rewrite rules in Mathematica's global rewrite table is similar to looking up the closure in a dictionary keyed on its function name. Such a dictionary can be constructed like this in F#:

> open System.Collections.Generic;; > let fns = Dictionary<string, (obj -> obj)>(dict["sqrt", unbox >> sqrt >> box]);; 

This is similar to the DownValues data structure in Mathematica, except that we aren't searching multiple resulting rules for the first to match on the function arguments.

The program then becomes:

> Array.map (fun x -> fns.["sqrt"] (box x)) xs;; Real: 00:00:00.044, CPU: 00:00:00.031, GC gen0: 0, gen1: 0, gen2: 0 ... 

Note that we get a similar 10× performance degradation due to the hash table lookup in the inner loop.

An alternative would be to store the DownValues associated with a symbol in the symbol itself in order to avoid the hash table lookup.

We can even write a complete term rewriter in just a few lines of code. Terms may be expressed as values of the following type:

> type expr = | Float of float | Symbol of string | Packed of float [] | Apply of expr * expr [];; 

Note that Packed implements Mathematica's packed lists, i.e. unboxed arrays.

The following init function constructs a List with n elements using the function f, returning a Packed if every return value was a Float or a more general Apply(Symbol "List", ...) otherwise:

> let init n f = let rec packed ys i = if i=n then Packed ys else match f i with | Float y -> ys.[i] <- y packed ys (i+1) | y -> Apply(Symbol "List", Array.init n (fun j -> if j<i then Float ys.[i] elif j=i then y else f j)) packed (Array.zeroCreate n) 0;; val init : int -> (int -> expr) -> expr 

The following rule function uses pattern matching to identify expressions that it can understand and replaces them with other expressions:

> let rec rule = function | Apply(Symbol "Sqrt", [|Float x|]) -> Float(sqrt x) | Apply(Symbol "Map", [|f; Packed xs|]) -> init xs.Length (fun i -> rule(Apply(f, [|Float xs.[i]|]))) | f -> f;; val rule : expr -> expr 

Note that the type of this function expr -> expr is characteristic of term rewriting: rewriting replaces expressions with other expressions rather than reducing them to values.

Our program can now be defined and executed by our custom term rewriter:

> rule (Apply(Symbol "Map", [|Symbol "Sqrt"; Packed xs|]));; Real: 00:00:00.049, CPU: 00:00:00.046, GC gen0: 24, gen1: 0, gen2: 0 

We've recovered the performance of Map[Sqrt, xs] in Mathematica!

We can even recover the performance of Sqrt[xs] by adding an appropriate rule:

| Apply(Symbol "Sqrt", [|Packed xs|]) -> Packed(Array.map sqrt xs) 

I wrote an article on term rewriting in F#.

1 Comment

Very interesting - +1. Actually, OCaml is on my wish list to learn (perhaps F# as well, but way lower in the list). One additional thing to take into account is that Map auto-compiles (JIT) a function on a list of length larger than SystemOptions["CompileOptions"->"MapCompileLength"], which is 100 by default. This is what really explains the 10-fold difference you noted.
6

Some measurements

Based on @gdelfino answer and comments by @rcollyer I made this small program:

j = # # + # # &; g[x_] := x x + x x ; h = Function[{x}, x x + x x ]; anon = Table[Timing[Do[ # # + # # &[i], {i, k}]][[1]], {k, 10^5, 10^6, 10^5}]; jj = Table[Timing[Do[ j[i], {i, k}]][[1]], {k, 10^5, 10^6, 10^5}]; gg = Table[Timing[Do[ g[i], {i, k}]][[1]], {k, 10^5, 10^6, 10^5}]; hh = Table[Timing[Do[ h[i], {i, k}]][[1]], {k, 10^5, 10^6, 10^5}]; ListLinePlot[ {anon, jj, gg, hh}, PlotStyle -> {Black, Red, Green, Blue}, PlotRange -> All] 

The results are, at least for me, very surprising: alt text

Any explanations? Please feel free to edit this answer (comments are a mess for long text)

Edit

Tested with the identity function f[x] = x to isolate the parsing from the actual evaluation. Results (same colors):

alt text

Note: results are very similar to this Plot for constant functions (f[x]:=1);

2 Comments

no edit capability. But, your anon function is not the same as g, h, and j, so of course it has a different result.
@rcoyller Tnx! I mixed up two versions. Corrected now. There is still the blue line (Function ...) to explain ...
4

Pattern matching seems faster:

In[1]:= g[x_] := x*x In[2]:= h = Function[{x}, x*x]; In[3]:= Do[h[RandomInteger[100]], {1000000}] // Timing Out[3]= {1.53927, Null} In[4]:= Do[g[RandomInteger[100]], {1000000}] // Timing Out[4]= {1.15919, Null} 

Pattern matching is also more flexible as it allows you to overload a definition:

In[5]:= g[x_] := x * x In[6]:= g[x_,y_] := x * y 

For simple functions you can compile to get the best performance:

In[7]:= k[x_] = Compile[{x}, x*x] In[8]:= Do[k[RandomInteger[100]], {100000}] // Timing Out[8]= {0.083517, Null} 

7 Comments

I edited my answer to show that maximum performance is achieved using a compiled function. But this is not always possible.
From The Help: * The number of times and the order in which objects are evaluated by Compile may be different from ordinary Mathematica code. * AND * Compiled code does not handle numerical precision and local variables in the same way as ordinary Mathematica code. *
@rcollyer In my machine j = # # &; is consistently faster than both
@belisarius, interesting. my machine also. So, it appears the use of the word Function introduces some overhead, but a pure function does much better. on my machine setting j=# #& and using j introduces some overhead v using # # directly.
@rcollyer: Coming very late to the party, but still... The reason why Function[x,...] is slower is that it is a lexical scoping construct that must care about possible variable name conflicts. The reason why pattern-matching is faster is that, while it is also a scoping construct, it does not care about the name conflicts in inner scoping constructs. The reason that #-& notation Function is the fastest is that it is not a lexical scoping construct in the sense that the others are, in other words it has the least number of checks to do.
|
3

You can use function recordSteps in previous answer to see what Mathematica actually does with Functions. It treats it just like any other Head. IE, suppose you have the following

f = Function[{x}, x + 2]; f[2] 

It first transforms f[2] into

Function[{x}, x + 2][2] 

At the next step, x+2 is transformed into 2+2. Essentially, "Function" evaluation behaves like an application of pattern matching rules, so it shouldn't be surprising that it's not faster.

You can think of everything in Mathematica as an expression, where evaluation is the process of rewriting parts of the expression in a predefined sequence, this applies to Function like to any other head

3 Comments

A late comment... The evaluation of Function may resemble pattern-matching, but is nevertheless fundamentally different. The reason for what you observed is that Function[{x},x+2] is indeed an expression, and that heads are evaluated before anything else. But it is not possible to explain how Function[x,Head[Unevaluated[x]],HoldAll][Print[1]] can ever return Print, if we consider Function as just a built-in symbol with SubValues. Function is independent from the pattern-matching mechanism, in this sense. Regarding the speed issue, see my comment for the other answer.
@Leonid Ah, there's no way to specify HoldAll for SubValues instead of DownValues, is there?
Yep, that's right. Could emulate, but not easy. When we use SetDelayed, the symbol is so-called symbolic head. For f[g,h][x,y] it is f. So, when f is HoldAll, it will hold g,h, but not x,y. So, if we assume that Function is just a head with built-in SubValues, there is no way to explain how it can hold the arguments passed to it. There is in fact a way to reproduce this feature with rules, using the evaluation Stack, but this is neither trivial nor robust, and surely not the way Function is implemented. It seems to be deeply wired in, one of the "magic symbols".
0

A thorough benchmark comparing all kinds of pattern matching with short names versus long names:

 164ns #1 & 167ns OddQ 394ns OddQ[#1] & 741ns If[OddQ[#1], 0, #1] & 1144ns Function[x, If[OddQ[x], 0, x]] 1139ns Function[electroencephalographically, If[OddQ[electroencephalographically], 0, electroencephalographically]] 815ns fn1 := If[OddQ[#1], 0, #1] & 1207ns fn2 := Function[x, If[OddQ[x], 0, x]] 1210ns fn3 := Function[electroencephalographically, If[OddQ[electroencephalographically], 0, electroencephalographically]] 864ns fn4[x_] := If[OddQ[x], 0, x] 844ns fn5[electroencephalographically_] := If[OddQ[electroencephalographically], 0, electroencephalographically] 829ns fn6[(x_)?OddQ] := 0 818ns fn7[(electroencephalographically_)?OddQ] := 0 836ns fn8[x_ /; OddQ[x]] := 0 839ns fn9[electroencephalographically_ /; OddQ[electroencephalographically]] := 0 845ns fnA[x_] /; OddQ[x] := 0 827ns fnB[electroencephalographically_] /; OddQ[electroencephalographically] := 0 1321ns fnC[x_] := 0 /; OddQ[x] 1333ns fnD[electroencephalographically_] := 0 /; OddQ[electroencephalographically] 1461ns fnE[(x_)?OddQ] := With[{v = x}, 0] 1485ns fnF[(electroencephalographically_)?OddQ] := With[{v = electroencephalographically}, 0] 2620ns fnG[x_] := With[{v = x}, 0 /; OddQ[v]] 2584ns fnH[electroencephalographically_] := With[{v = electroencephalographically}, 0 /; OddQ[v]] 

Code tested with Mathematica 14.3.0 on Linux x86_64:

ClearAll["`*"];ClearSystemCache[]; Attributes@tstF={HoldAll,Listable};Attributes@strF={HoldAllComplete}; valF[st_]:=ToString[HoldForm@st/.{HoldPattern->HoldForm,RuleDelayed->SetDelayed},TextForm] strF[ex_:Except@String_]:=valF@FirstCase[Quiet@{DownValues@ex,OwnValues@ex},_RuleDelayed,HoldForm@ex,2] (tstF[fn_]:=StringPadLeft[ToString@Round@Median@Table[300030First@Timing@CompoundExpression@##,31]<>"ns",6]<>" "<>strF@fn)&@@(fn/@Range@3333); fn1=If[OddQ[#],0,#]&; fn2=x|->If[OddQ[x],0,x]; fn3=electroencephalographically|->If[OddQ[electroencephalographically],0,electroencephalographically]; fn4[x_]:=If[OddQ[x],0,x] fn5[electroencephalographically_]:=If[OddQ[electroencephalographically],0,electroencephalographically] fn6[x_?OddQ]:=0;fn6[x_?EvenQ]:=x fn7[electroencephalographically_?OddQ]:=0;fn7[electroencephalographically_?EvenQ]:=electroencephalographically fn8[x_/;OddQ@x]:=0;fn8[x_/;EvenQ@x]:=x fn9[electroencephalographically_/;OddQ@electroencephalographically]:=0;fn9[electroencephalographically_/;EvenQ@electroencephalographically]:=electroencephalographically fnA[x_]/;OddQ@x:=0;fnA[x_]/;EvenQ@x:=x fnB[electroencephalographically_]/;OddQ@electroencephalographically:=0;fnB[electroencephalographically_]/;EvenQ@electroencephalographically:=electroencephalographically fnC[x_]:=0/;OddQ@x;fnC[x_]:=x/;EvenQ@x fnD[electroencephalographically_]:=0/;OddQ@electroencephalographically;fnD[electroencephalographically_]:=electroencephalographically/;EvenQ@electroencephalographically fnE[x_?OddQ]:=With[{v=x},0];fnE[x_?EvenQ]:=With[{v=x},v] fnF[electroencephalographically_?OddQ]:=With[{v=electroencephalographically},0];fnF[electroencephalographically_?EvenQ]:=With[{v=electroencephalographically},v] fnG[x_]:=With[{v=x},0/;OddQ@v];fnG[x_]:=With[{v=x},v/;EvenQ@x] fnH[electroencephalographically_]:=With[{v=electroencephalographically},0/;OddQ@v];fnH[electroencephalographically_]:=With[{v=electroencephalographically},v/;EvenQ@electroencephalographically] RawBoxes@TextData@ToString@Column@tstF@{#&,OddQ,OddQ@#&,Evaluate@fn1,Evaluate@fn2,Evaluate@fn3,fn1,fn2,fn3,fn4,fn5,fn6,fn7,fn8,fn9,fnA,fnB,fnC,fnD,fnE,fnF,fnG,fnH} 

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.