一个高效的向量幂运算方法

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 比其他方法快两倍。然而，当我计算 x.^6 时，(x.*x.*x).^2 和 x.*x.*x.*x.*x.*x 的差异非常小（而其他选项则较慢）。

有没有办法确定对向量进行幂运算的最有效方式？为什么性能差别如此之大？

- user2041376

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的次幂进行一次幂运算。这个技巧适用于所有Legendre多项式（在奇次多项式中，x2被替换为x）。

- Luis Mendo

Answer 2

这里有一些想法： power(x,4) 和 x.^4 是等价的（只需阅读文档）。 x.*x.*x.*x 可能被优化为类似于 x.^2.^2 的东西。

---

x.^2.^2 可能会被解释为：对每个元素进行平方（快速），再次对结果进行平方（再次快速）。

x.^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倍的差异。但请注意，以秒为单位的时间差仍然很小，因此对于大多数实际应用程序，我只会选择简单的语法。

- Dennis Jaheruddin

---

对于 x.*x.*x.*x 可能进行的优化表现出奇怪的行为。我尝试了 x.*.x.* ... .*x，其中 "x" 的数量从 2 到 8 不等，时间增加得多少不一。我本来期望会有波动；例如 "8" 情况（=> x.^2.^2.^2：三个幂运算）应该比 "7"（=> 更多的幂运算）花费更少的时间。 - Luis Mendo

@LuisMendo 我不知道如何验证，但我可以想象它只执行1步（没有嵌套优化）。对于7，它将简化为类似于：x.^2*x.^2*x.^2.*x，这不会比8的x.^2*x.^2*x.^2.*x.^2慢。如果用这种方式做8比做7更快，Mathworks可能已经在幂函数中实现了这种优化。 - Dennis Jaheruddin

是的，那可能就是解释：没有嵌套。 - Luis Mendo

@DennisJaheruddin，我认为你是正确的。看看我的答案（当你回答时我正在撰写）- 嵌套对于16次方来说出人意料的慢2倍。 - mbauman

Answer 3

似乎 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

- mbauman