获取向量幂的有效方法

Question

我编写了一个代码，它在数值上使用了最高 n 阶的勒让德多项式。例如：

....
case 8 
p = (6435*x.^8-12012*x.^6+6930*x.^4-1260*x.^2+35)/128; return
case 9 
...

如果向量x 很长，这可能会变慢。我看到x.^4 和x.*x.*x.*x 之间存在性能差异，并认为我可以使用它来改进我的代码。我使用了timeit 并发现：

x=linspace(0,10,1e6);
f1= @() power(x,4)
f2= @() x.4;
f3= @() x.^2.^2
f4= @() x.*x.*x.*x

f4 比其他人快 因素 2。但是，当我转到 x.^6 时，(x.*x.*x).^2 和 x.*x.*x.*x.*x.*x 之间几乎没有区别（而所有其他选项都比较慢）。

有没有办法告诉我们什么是获取向量幂的最有效方法？ 你能解释一下为什么会有这么大的性能差异吗？

Answer 1

这并不完全是您问题的答案，但它可能会解决您的问题：

x2 = x.*x; % or x.^2 or power(x,2), whichever is most efficient
p = ((((6435*x2-12012)*x2+6930)*x2-1260)*x2+35)/128

这样你只做一次幂，而且只用指数 2。这个技巧可以应用于所有勒让德多项式（在奇数次多项式中，一个 x2 被 x 替换）。

Answer 2

似乎 Mathworks 在其幂函数中有特殊的大小写正方形（不幸的是，它都是我们看不到的内置封闭源代码）。在我对 R2013b 的测试中，.^、power 和 realpow 似乎使用相同的算法。对于正方形，我相信他们已将其特例为x.*x。

1.0x (4.4ms):   @()x.^2
1.0x (4.4ms):   @()power(x,2)
1.0x (4.5ms):   @()x.*x
1.0x (4.5ms):   @()realpow(x,2)
6.1x (27.1ms):  @()exp(2*log(x))

对于立方体，情况有所不同。它们不再是特殊情况。同样，.^、power 和 realpow 都相似，但这次要慢得多：

1.0x (4.5ms):   @()x.*x.*x
1.0x (4.6ms):   @()x.*x.^2
5.9x (26.9ms):  @()exp(3*log(x))
13.8x (62.3ms): @()power(x,3)
14.0x (63.2ms): @()x.^3
14.1x (63.7ms): @()realpow(x,3)

让我们跳到 16 次方，看看这些算法如何扩展：

1.0x (8.1ms):   @()x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x
2.2x (17.4ms):  @()x.^2.^2.^2.^2
3.5x (27.9ms):  @()exp(16*log(x))
7.9x (63.8ms):  @()power(x,16)
7.9x (63.9ms):  @()realpow(x,16)
8.3x (66.9ms):  @()x.^16

所以：.^、power 和 realpow 都在关于指数的恒定时间内运行，除非它是特殊情况（-1 也似乎是特殊情况）。使用exp(n*log(x)) 技巧也是关于指数的恒定时间，而且速度更快。唯一的结果我不太明白为什么重复平方比乘法慢。

正如预期的那样，将 x 的大小增加 100 倍会增加所有算法的时间。

那么，这个故事的寓意是什么？使用标量整数指数时，请始终自己进行乘法运算。 power 和朋友们有很多聪明之处（指数可以是浮点数、向量等）。唯一的例外是 Mathworks 为您完成优化的地方。在 2013b 中，似乎是x^2 和x^(-1)。希望随着时间的推移他们会增加更多。但是，一般来说，求幂很困难，而乘法很容易。在性能敏感的代码中，我认为总是输入x.*x.*x.*x 不会出错。 （当然，在你的情况下，请遵循 Luis 的建议，并利用每个学期的中间结果！）

function powerTest(x)

f{1} = @() x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x;
f{2} = @() x.^2.^2.^2.^2;
f{3} = @() exp(16.*log(x));
f{4} = @() x.^16;
f{5} = @() power(x,16);
f{6} = @() realpow(x,16);

for i = 1:length(f)
    t(i) = timeit(f{i});
end

[t,idxs] = sort(t);
fcns = f(idxs);

for i = 1:length(fcns)
    fprintf('%.1fx (%.1fms):\t%s\n',t(i)/t(1),t(i)*1e3,func2str(fcns{i}));
end

Answer 3

以下是一些想法：

power(x,4) 和 x.^4 是等效的（只需阅读文档）。

x.*x.*x.*x 可能已优化为类似x.^2.^2

---

x.^2.^2 可能被评估为：取每个元素的平方（快速），然后再次取其平方（再次快速）。

x.^4 可能直接评估为：取每个元素的四次方（慢）。

看到 2 次快速操作比 1 次慢速操作花费的时间更少，这并不奇怪。太糟糕了，在幂 4 的情况下没有执行优化，但也许它并不总是有效或需要付出代价（输入检查、内存？）。

---

关于时间安排：实际上差异远不止 2 倍！

当您现在在函数中调用它们时，在每种情况下都会增加函数开销，从而使相对差异更小：

y=x;tic,power(x,4);toc
y=x;tic,x.^4;toc
y=x;tic,x.^2.^2;toc
y=x;tic,x.*x.*x.*x;toc

将给予：

Elapsed time is 0.034826 seconds.
Elapsed time is 0.029186 seconds.
Elapsed time is 0.003891 seconds.
Elapsed time is 0.003840 seconds.

因此，这几乎是 10 倍的差异。但是，请注意，以秒为单位的时间差异仍然很小，因此对于大多数实际应用程序，我只会使用简单的语法。