Mathematica中函数与模式匹配的性能差异

我对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 . 我建议你看看 introduction 利昂尼德希夫林的优秀的书,其中这些概念解释了一些细节。

我有时会建造一个 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公司 所以搜索全局重写规则是非常昂贵的,而术语重写是昂贵的。

如果是,为什么 Sqrt公司 快吗?您可能会期望2倍而不是10倍的速度减慢,因为我们已经替换了一个查找 Sqrt公司 两次查找 sqrt公司 和 Sqrt公司 但这并不是因为 Sqrt公司 具有内置函数的特殊状态,该函数将与Mathematica术语重写器本身的核心相匹配,而不是通过通用全局重写表。

其他人描述过相似功能之间的性能差异要小得多。我相信这些情况下的性能差异只是Mathematica内部的精确实现中的微小差异。Mathematica最大的问题是全局重写表。特别是,这是Mathematica与传统术语级解释器的不同之处。

通过编写mini 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。但这是意料之中的,因为这段代码已经被.NET编译成机器代码,以便快速计算。

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

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

这与 DownValues Mathematica中的数据结构,除了我们没有搜索多个结果规则以使第一个规则与函数参数匹配之外。

然后程序变为:

> 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倍。

另一种方法是存储 下降值 与符号本身中的符号关联,以避免哈希表查找。

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

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

请注意 Packed 实现Mathematica的压缩列表,即未绑定数组。

以下 init 函数构造 List 具有 n 使用函数的元素 f ,返回 拥挤的 如果每个返回值都是 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<i then Float ys.[i]
              elif j=i then y
              else f j))
    packed (Array.zeroCreate n) 0;;
val init : int -> (int -> expr) -> expr

以下 rule 函数使用模式匹配来标识它可以理解的表达式,并用其他表达式替换它们:

> 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

注意这个函数的类型 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

我们已经恢复了 Map[Sqrt, xs] 在Mathematica!

我们甚至可以恢复 Sqrt[xs] 通过添加适当的规则:

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

我写了一篇关于 term rewriting in F# .

一些测量

根据@gdelfino的回答和@rcollyer的评论,我制作了这个小程序:

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]

结果,至少对我来说,非常令人惊讶:

有什么解释吗?请随意编辑此答案(对于长文本,注释是混乱的)

编辑

使用标识函数f[x]=x进行测试,以将解析与实际计算隔离开来。结果(相同颜色):

注:对于常数函数(f[x]:=1),结果与此图非常相似;

模式匹配似乎更快:

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}

模式匹配也更灵活,因为它允许您重载定义:

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

对于简单的函数,可以编译以获得最佳性能:

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

您可以使用前面的函数recordSteps answer 看看Mathematica对函数到底做了什么。它对待它就像对待其他人一样。也就是说,假设你有

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

它首先将f[2]转换为

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

下一步, x+2 转化为 2+2 . 从本质上说,“函数”求值的行为类似于模式匹配规则的应用程序,因此它不快也就不足为奇了。

您可以将Mathematica中的所有内容都看作一个表达式,其中求值是在 predefined sequence ,这适用于类似于任何其他头部的功能