解决不用for循环的广播错误，加速代码

Question

我可能误解了 Python 中广播的工作原理，但我仍然运行出错。

scipy 提供了一些接受两个参数的 "special functions"，特别是 eval_XX(n, x[,out]) 函数。参见 http://docs.scipy.org/doc/scipy/reference/special.html

我的程序使用了许多正交多项式，所以我必须在不同的点计算这些多项式。举个具体的例子scipy.special.eval_hermite(n, x, out=None)。

我希望 x 参数是矩阵形状 (50, 50)。然后，我想在多个点上评估该矩阵的每个条目。让我们将 n 定义为一个 numpy 数组 narr = np.arange(10)（我们将 numpy 导入为 np，即 import numpy as np）。

所以，调用

scipy.special.eval_hermite(narr, matrix)

应该return厄米多项式H_0(matrix), H_1(matrix), H_2(matrix)，等等。每个H_X(matrix)的形状是(50,50)，原始输入矩阵的形状。

然后，我想对这些值求和。所以，我打电话给

matrix1 = np.sum( [scipy.eval_hermite(narr, matrix)], axis=0 )

但是我收到一个广播错误！

ValueError: operands could not be broadcast together with shapes (10,) (50,50)

我可以用 for 循环解决这个问题，即

matrix2 = np.sum( [scipy.eval_hermite(i, matrix) for i in narr], axis=0)

这给了我正确的答案，输出 matrix2.shape = (50,50)。但是使用这个 for 循环会大大降低我的代码速度。请记住，我们正在处理矩阵的条目。

有没有不用 for 循环的方法？

Answer 1

这些函数的文档很简陋，而且很多代码都是经过编译的，所以这只是基于实验：

special.eval_hermite(n, x, out=None)

n 显然是一个标量或整数数组。 x 可以是浮点数数组。

special.eval_hermite(np.ones(5,int)[:,None],np.ones(6)) 给了我 (5,6) 结果。这与我从 np.ones(5,int)[:,None] * np.ones(6).

得到的形状相同

np.ones(5,int)[:,None] 是一个 (5,1) 数组，np.ones(6) 是一个 (6,)，就此而言，它等同于 (1,6)。两者都可以展开为(5,6).

据我所知，这些 special 函数中的广播规则与 *.

等运算符相同

由于 special.eval_hermite(nar[:,None,None], x) 产生 (10,50,50)，您只需将 sum 应用于轴 0 即可产生 (50,50).

special.eval_hermite(nar[:,Nar,Nar], x).sum(axis=0)

就像我之前写的那样，相同的广播（和求和）规则适用于此 hermite 就像它们适用于 *.

等基本操作一样

Answer 2

eval_hermite 用 x 广播 n，然后在每个点计算 H_n(x)。因此，输出形状将是广播 n 和 x 的结果。所以，如果你想完成这项工作，你必须让 n 和 x 具有兼容的形状：

import scipy.special as ss
import numpy as np
matrix = np.ones([100,100]) # example
narr = np.arange(10) # example
ss.eval_hermite(narr[:,None,None], matrix).shape # => (10, 100, 100)

但请注意，这实际上可能会更快：

out = np.zeros_like(matrix)
for n in narr:
    out += ss.eval_hermite(n, matrix)

在测试中，它似乎比上面的 np.sum(...) 快 5-10%。

解决不用for循环的广播错误，加速代码

Solve broadcasting error without for loop, speed up code

python

for-loop

numpy

matrix

scipy