系数函数速度较慢

Question

请考虑以下内容：

Clear[x]
expr = Sum[x^i, {i, 15}]^30;

CoefficientList[expr, x]; // Timing
Coefficient[Expand@expr, x, 234]; // Timing
Coefficient[expr, x, 234]; // Timing

{0.047，空值}

{0.047，空值}

{4.93，空值}

帮助状态：

Coefficient 可以处理 expr 是否以展开形式明确给出的情况。

为什么在最后一种情况下 Coefficient 需要如此缓慢？

- Mr.Wizard

---

也许 Coefficient 算法会以速度为代价来换取能够处理极长展开式的表达式的能力？顺便说一下，你的电脑比我的快4.5倍。 - Szabolcs

@Szabolcs 我想这是有道理的；我会尝试测试一下。如果是这样的话，我希望能有一个更智能的方法选择过程。如果您正在使用版本8，您也可以尝试在版本7中进行测试吗？我怀疑至少有些东西在8中会变慢。 - Mr.Wizard

@Szabolcs 我可以证实你的假设，至少在 (1 + x)^50000 上，Coefficient 更加节省内存。有没有什么合理的方法可以让一个调用 Coefficient 的通用函数更快一些？是否有某种半展开形式或 Method 选项可以在这些选项之间给我一个平衡点？ - Mr.Wizard

表达式越大情况就越糟。 当expr = Sum[x ^ i，{i，30}] ^ 20;时，它们之间有670的因素。 就内存使用而言，我没有注意到太大的区别，但系统监视器不是非常敏感。 - Sjoerd C. de Vries

你看到另一个奇怪的版本差异了吗（在8.0.x和8.0.4之间）这里？Belisarius的解决方案在8.0.4上无法正常工作，显然是由于Radon[]的更改。 - Szabolcs

Answer 1

这里有一个可能让你的代码更快的技巧，但我不能保证它总是能正确运行：

ClearAll[withFastCoefficient];
SetAttributes[withFastCoefficient, HoldFirst];
withFastCoefficient[code_] :=
   Block[{Binomial},
     Binomial[x_, y_] := 10 /; ! FreeQ[Stack[_][[-6]], Coefficient];
     code]

使用它，我们得到：

In[58]:= withFastCoefficient[Coefficient[expr,x,234]]//Timing
Out[58]= {0.172,3116518719381876183528738595379210}

这个想法是，Coefficient 在内部使用 Binomial 估计项数，如果项数小于 1000，则扩展（调用 Expand），您可以通过使用 Trace[..., TraceInternal->True] 进行检查。当它不扩展时，它会计算许多大量由零主导的系数列表的求和，这显然比扩展慢，适用于一系列表达式。我的做法是欺骗 Binomial 返回一个小数（10），但我还尝试使其仅影响 Coefficient 内部调用的 Binomial：

In[67]:= withFastCoefficient[f[Binomial[7,4]]Coefficient[expr,x,234]]
Out[67]= 3116518719381876183528738595379210 f[35] 

我不能保证在代码的其他地方计算Binomial时不会出现错误。

编辑

当然，一种更安全的替代方法是使用Villegas - Gayley技巧重新定义Coefficient，展开其中的表达式并再次调用：

Unprotect[Coefficient];
Module[{inCoefficient},
   Coefficient[expr_, args__] :=
      Block[{inCoefficient = True},
         Coefficient[Expand[expr], args]] /; ! TrueQ[inCoefficient]
];
Protect[Coefficient];

编辑2

我的第一个建议有一个优点，即我们定义了一个宏，可以在本地修改函数的属性，但缺点是不安全。我的第二个建议更加安全，但会在全局范围内修改Coefficient，因此它将一直扩展，直到我们删除该定义。我们可以借助Internal`InheritedBlock来实现两全其美，它可以创建给定函数的本地副本。以下是代码：

ClearAll[withExpandingCoefficient];
SetAttributes[withExpandingCoefficient, HoldFirst];
withExpandingCoefficient[code_] :=
   Module[{inCoefficient},
      Internal`InheritedBlock[{Coefficient},
         Unprotect[Coefficient];
         Coefficient[expr_, args__] :=
            Block[{inCoefficient = True},
               Coefficient[Expand[expr], args]] /; ! TrueQ[inCoefficient];
         Protect[Coefficient];
         code
      ]
   ];

使用方式与第一种情况类似：

In[92]:= withExpandingCoefficient[Coefficient[expr,x,234]]//Timing
Out[92]= {0.156,3116518719381876183528738595379210}

然而，主要的Coefficient函数保持不变：

In[93]:= DownValues[Coefficient]
Out[93]= {}

- Leonid Shifrin

---

对于逆向工程和分析的工作，我给予+1的认可，但这难道不是与手动展开相同，只是不安全吗？或者它仍然比手动展开更节省内存（即是否对整个表达式进行扩展或不扩展的内部决策，还是仅针对部分）？ - Szabolcs

@Szabolcs 我刚刚编辑了一下，在看到你的评论之前添加了另一种可能性（扩展）。从内存方面来说，它可能是相同的，而且当系数展开时，它会展开所有表达式（据我所知）。我的第一个解决方案更多是为了说明。我很快就会添加基于Gayley-Villegas技巧的第二个解决方案。 - Leonid Shifrin

Answer 2

Coefficient 只有在认为有必要扩展时才会进行扩展，这确实避免了内存爆炸。我相信自从3版本以来就一直是这样的（我记得我大约在1995年左右开始使用它）。

避免扩展也可能更快。下面是一个简单的例子。

In[28]:= expr = Sum[x^i + y^j + z^k, {i, 15}, {j, 10}, {k, 20}]^20;

In[29]:= Coefficient[expr, x, 234]; // Timing

Out[29]= {0.81, Null}

但是这个在第8版中似乎会卡住，而且在开发版的Mathematica中至少需要半分钟（因为Expand已经变了）。

Coefficient[Expand[expr], x, 234]; // Timing

可能需要添加一些启发式方法来查找不会爆炸的单变量。但似乎不是非常重要的事情。

Daniel Lichtblau

- Daniel Lichtblau

Answer 3

expr = Sum[x^i, {i, 15}]^30; 

scoeff[ex_, var_, n_] /; PolynomialQ[ex, var] := 
   ex + O[var]^(n + 1) /. 
    Verbatim[SeriesData][_, 0, c_List, nmin_, nmax_, 1] :> 
     If[nmax - nmin != Length[c], 0, c[[-1]]]; 

Timing[scoeff[expr, x, 234]]

看起来速度相当快。

- Rolf Mertig

---

+1. 对于较大的表达式，它实际上比基于展开表达式的解决方案快几倍。此外，对于较小指数系数，它的速度显著更快，因为它似乎只扩展到指定项（幂）的表达式部分。 - Leonid Shifrin

SeriesData 在 Mathematica 中的实现非常出色。使用 SeriesData 可以轻松地完成任何任务，且速度极快（甚至比 FORM（http://www.nikhef.nl/~form/）还要快）。 - Rolf Mertig

Rolf Mertig，我被迫接受这个答案，因为它简洁而有效。 但是，我并不真正理解它。您能以可能最简单的方式解释一下吗？ - Mr.Wizard

具体来说，我不明白的是：ex + O[var]^(n + 1)。 - Mr.Wizard

这只是一个泰勒展开，在这里有意义，因为x的高次幂不需要。加速来自于SeriesData在Mathematica中高效实现的事实。另请参阅：http://reference.wolfram.com/mathematica/ref/O.html - Rolf Mertig

谢谢。对不起，我取消了您的答案，但我认为我找到了一种更简单的方法，受到您的启发。我会等待同行对该答案的审查。 - Mr.Wizard

Answer 4

在跟随Rolf Mertig的答案做出一些实验后，这似乎是处理像Sum[x^i, {i, 15}]^30这样类型表达式的最快方法：

SeriesCoefficient[expr, {x, 0, 234}]

- Mr.Wizard

---

请注意，这个回答与您最初提出的问题略有不同：您问的是Coefficient缓慢的原因（重点在于Coefficient和“为什么缓慢”），而不是可能更快的替代方法。 - Leonid Shifrin

@Leonid 非常好的观点。在阅读了各种答案后，我忘记了这一点。我想Daniel的回答最直接;我应该接受吗？此外，你认为这种方法怎么样？它有任何明显的问题吗？ - Mr.Wizard

抱歉，我无法建议您接受哪个答案 - 这完全取决于您。丹尼尔的回答非常好，简洁明了。至于方法 - 我没有看到明显的问题，但是很遗憾现在没有时间去尝试它。 - Leonid Shifrin

@Leonid，你听起来很圆滑。好的。 - Mr.Wizard