是否有“增强的”numpy/scipy dot 方法？

Posted 2023-03-12

【中文标题】是否有“增强的”numpy/scipy dot 方法？

【英文标题】：Is there an "enhanced" numpy/scipy dot method?

【发布时间】：2012-03-17 17:44:02

【问题描述】：

问题

我想使用 numpy 或 scipy 计算以下内容：

Y = A**T * Q * A

其中A 是m x n 矩阵，A**T 是A 的转置，Q 是m x m 对角矩阵。

由于Q 是对角矩阵，我只将其对角元素存储为向量。

Y 的求解方法

目前我能想到两种计算Y的方法：

Y = np.dot(np.dot(A.T, np.diag(Q)), A) 和 Y = np.dot(A.T * Q, A)。

显然选项 2 比选项 1 更好，因为不必使用 diag(Q) 创建真正的矩阵（如果这是 numpy 真正做的......） 然而，这两种方法都有一个缺陷，即必须分配比实际需要更多的内存，因为A.T * Q 和np.dot(A.T, np.diag(Q)) 必须与A 一起存储才能计算Y。

问题

在 numpy/scipy 中是否有一种方法可以消除不必要的额外内存分配，在这种情况下您只会传递两个矩阵 A 和 B（在我的情况下 B 是 A.T）和一个加权向量Q 一起来吗？

【参考方案1】：

（w/r/t OP 的最后一句话：我不知道这种 numpy/scipy 方法，但是 w/r/t OP 标题中的问题（即改进NumPy 点性能）下面的内容应该会有所帮助。换句话说，我的回答是针对提高大多数步骤的性能，这些步骤包括你的 Y) 函数。

首先，这应该会给您带来比普通 NumPy dot 方法明显的提升：

>>> from scipy.linalg import blas as FB
>>> vx = FB.dgemm(alpha=1., a=v1, b=v2, trans_b=True)

请注意，v1、v2 这两个数组在 C_FORTRAN 顺序中都

您可以通过数组的 flags 属性访问 NumPy 数组的字节顺序，如下所示：

>>> c = NP.ones((4, 3))
>>> c.flags
      C_CONTIGUOUS : True          # refers to C-contiguous order
      F_CONTIGUOUS : False         # fortran-contiguous
      OWNDATA : True
      MASKNA : False
      OWNMASKNA : False
      WRITEABLE : True
      ALIGNED : True
      UPDATEIFCOPY : False

要更改其中一个数组的顺序，使两者对齐，只需调用 NumPy 数组构造函数，传入数组并将适当的 order 标志设置为 True

>>> c = NP.array(c, order="F")

>>> c.flags
      C_CONTIGUOUS : False
      F_CONTIGUOUS : True
      OWNDATA : True
      MASKNA : False
      OWNMASKNA : False
      WRITEABLE : True
      ALIGNED : True
      UPDATEIFCOPY : False

您可以通过利用数组顺序对齐来进一步优化减少因复制原始数组而导致的过多内存消耗。

但是为什么数组在传递给dot之前要被复制呢？

点积依赖于 BLAS 运算。这些操作需要以 C 连续顺序存储的数组——正是这种约束导致数组被复制。

另一方面，转置确实不影响副本，但不幸的是以Fortran顺序返回结果：

因此，要消除性能瓶颈，您需要消除谓词数组复制步骤；要做到这一点，只需将两个数组都以 C 连续顺序传递给 dot*。

所以要计算 dot(A.T., A) 而不制作额外的副本：

>>> import scipy.linalg.blas as FB
>>> vx = FB.dgemm(alpha=1.0, a=A.T, b=A.T, trans_b=True)

总之，上面的表达式（连同谓词导入语句）可以代替点，提供相同的功能但性能更好

你可以像这样将该表达式绑定到一个函数：

>>> super_dot = lambda v, w: FB.dgemm(alpha=1., a=v.T, b=w.T, trans_b=True)

【讨论】：

这确实是访问 BLAS 例程的好方法。我相信我将来可以很好地利用它。但是，仍然有一个Q 要插入这里的某个地方...... :)

当然，我应该在上面的回答前加上“我不知道您在 OP 的 最后 句中描述的这种 numpy/scipy 方法，但是这是如何提高构成 Y 函数的大多数步骤的性能；如果这具有误导性，我们深表歉意（我的答案已相应编辑）。

@doug 速度呢？在我的测试中，无论我的矩阵是什么顺序，dgemm 似乎都比 np.dot 慢得多。 dgemm 应该比 numpy.dot 快吗？

你能解释一下'trans_b=True'是什么意思吗？参考文献没有对此的详细描述。

【参考方案2】：

我只是想把它放在 SO 上，但是这个拉取请求应该会有所帮助，并且不需要为 numpy.dot 提供单独的函数 https://github.com/numpy/numpy/pull/2730 这应该在 numpy 1.7 中可用

与此同时，我使用上面的示例编写了一个可以替换 numpy dot 的函数，无论数组的顺序如何，并正确调用 fblas.dgemm。 http://pastebin.com/M8TfbURi

希望这会有所帮助，

【讨论】：

您是否有任何示例说明 np.dot（在 numpy 1.8 中）何时/不复制 arg ？或者，我怎么知道？

【参考方案3】：

numpy.einsum 是您正在寻找的：

numpy.einsum('ij, i, ik -> jk', A, Q, A)

这不需要任何额外的内存（尽管通常 einsum 的工作速度比 BLAS 操作慢）