Mathematica中函数和模式匹配之间的性能差异

所以Mathematica与lisp的其他方言不同,因为它模糊了函数和宏之间的界限.在Mathematica中,如果用户想要编写数学函数,他们可能会使用模式匹配,f[x_]:= x*x而不是f=Function[{x},x*x]两者都会在调用时返回相同的结果f[x].我的理解是第一种方法相当于一个lisp宏,并且由于语法更简洁,我的经验受到青睐.

所以我有两个问题,执行函数与模式匹配/宏方法之间是否存在性能差异？虽然函数实际上转换为某些版本的宏以允许实现的功能,但我的一部分不会感到惊讶Listable.

我关心这个问题的原因是因为最近关于试图在大型程序中捕获Mathematica错误的问题(1) (2).如果大多数计算是根据函数定义的,那么在我看来,跟踪评估的顺序和错误发生的位置比在连续应用宏来重写输入后尝试捕获错误更容易/模式.

我理解Mathematica的方式是它是一个巨大的搜索替换引擎.所有函数,变量和其他赋值都基本上存储为规则,在评估期间,Mathematica会遍历此全局规则库并应用它们,直到结果表达式停止更改.

因此,您必须通过规则列表的次数越少,评估的速度就越快.看看使用什么Trace(使用gdelfino的函数g和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}

很明显为什么匿名函数是最快的,为什么使用Function引入额外的开销而不是简单SetDelayed.我建议在看引进列昂尼德·希夫林的优秀图书,其中这些概念进行了详细解释.

我偶尔会构建一个Dispatch包含我需要的所有函数的表,并将其手动应用到我的起始表达式中.与正常评估相比,这提供了显着的速度提升,因为Mathematica的内置功能都不需要与我的表达相匹配.

我的理解是第一种方法相当于一个lisp宏,并且由于语法更简洁,我的经验受到青睐.

并不是的.Mathematica是一个术语重写器,Lisp宏也是如此.

所以我有两个问题,执行函数与模式匹配/宏方法之间是否存在性能差异？

是.请注意,您从未真正在Mathematica中"执行函数".您只是应用重写规则将一个表达式更改为另一个表达式.

考虑将Sqrt函数映射到浮点数的压缩数组上.Mathematica中最快的解决方案是将Sqrt函数直接应用于打包数组,因为它恰好实现了我们想要的并且针对这种特殊情况进行了优化:

In[1] := N@Range[100000];

In[2] := Sqrt[xs]; // AbsoluteTiming

Out[2] = {0.0060000, Null}

我们可能会定义一个全局重写规则,其中包含sqrt[x]重写形式的术语,Sqrt[x]以便计算平方根:

In[3] := Clear[sqrt];
         sqrt[x_] := Sqrt[x];
         Map[sqrt, xs]; // AbsoluteTiming

Out[3] = {0.4800007, Null}

请注意,这比之前的解决方案慢约100倍.

或者,我们可以定义一个全局重写规则,sqrt用一个调用的lambda函数替换该符号Sqrt:

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

Out[4] = {0.0500000, Null}

请注意,这比之前的解决方案快约10倍.

为什么？因为慢速第二个解决方案是sqrt[x_] :> Sqrt[x]在内部循环中查找重写规则(对于数组的每个元素),而快速第三个解决方案查找Function[...]符号的值sqrt一次,然后重复应用该lambda函数.相比之下,最快的第一个解决方案是sqrt用C语言编写的循环调用.因此,搜索全局重写规则非常昂贵,并且术语重写很昂贵.

如果是这样,为什么要Sqrt快？你可能会想到一个2×放缓,而不是10×,因为我们已经更换了一个查找了Sqrt两个用于查找sqrt和Sqrt在内部循环,但事实并非如此,因为Sqrt有被内置的功能将在匹配的特殊地位Mathematica术语的核心是重写器本身,而不是通过通用的全局重写表.

其他人已经描述了类似功能之间更小的性能差异.我相信这些情况下的性能差异只是Mathematica内部的确切实现方面的微小差异.Mathematica最大的问题是全局重写表.特别是,这就是Mathematica与传统的术语级解释器不同的地方.

通过编写迷你Mathematica实现,您可以了解Mathematica的性能.在这种情况下,可以将上述解决方案编译为(例如)F#.可以像这样创建数组:

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

内置sqrt函数可以转换为闭包,并赋予map函数,如下所示:

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

这需要6毫秒,就像Sqrt[xs]在Mathematica中一样.但这是可以预料到的,因为这段代码已被JIT编译成.NET的机器代码,以便快速评估.

在Mathematica的全局重写表中查找重写规则类似于在键入其函数名的字典中查找闭包.这样的字典可以在F#中这样构造:

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

这类似于DownValuesMathematica中的数据结构,除了我们没有搜索多个结果规则以便第一个匹配函数参数.

该程序然后变为:

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

请注意,由于内部循环中的哈希表查找,我们得到了类似的10倍性能下降.

另一种方法是DownValues在符号本身中存储与符号相关联的符号,以避免哈希表查找.

我们甚至可以用几行代码编写一个完整的术语重写器.条款可以表示为以下类型的值:

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

请注意,Packed实现Mathematica的打包列表,即未装箱的数组.

以下init函数使用函数构造a Listwith n元素f,返回Packedif if每个返回值是a Float或更通用的Apply(Symbol "List", ...)否则:

> 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 (int -> expr) -> expr

> 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

> 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

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

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]