问numpy累积点积
EN

Stack Overflow用户

提问于 2019-12-26 11:09:19

回答 1查看 642关注 0票数 1

我有一个ndarray A，填充N个平方DxD矩阵(shape (N，D，D))。我想把它转换成一个形状相同的ndarray B，其中B=A和每一个i>0，Bi = np.dot(Bi-1，Ai)。虽然基本实现是显而易见的，但我想知道这个操作是否比for循环有更快的实现。

例如，让我描述另一种执行计算的方法：

B0.N/2=计算A.AN/2-1基本方式
BN/2.N=计算AN/2.a
返回np.concatenate((B0.N/2-1，np.dot(BN/2 .N))

重点是1和2可以并行完成，3是向量化操作，这种拆分可以根据需要对数组的每一半进一步应用。这让我怀疑是否存在比基本for循环更好的选项(例如，我建议的是实现/实际改进，还是另一个选项更好)。

非常感谢，

伊夫塔赫

编辑:用于大多数基本实现的代码，用于基准测试：

import numpy as np

def cumdot(A):
    B = np.empty(A.shape)
    B[0] = A[0]
    for i in range(1, A.shape[0]):  
        B[i] = B[i - 1] @ A[i]
    return B

Edit2:在numpy中，所有运行函数似乎都支持.accumulate() (这正是我所要做的)，matmul (它的行为就像一个点积)是一个通用的ufunc。这意味着matmul不是从两个标量到一个矩阵的函数，而是从两个矩阵到一个矩阵的函数，因此当函数累积存在时，调用matmul将引发一个异常，说明累积不能在具有签名的ufuncs上调用。如果这件事能成功的话，我也很想知道。

python

numpy

matrix

回答 1

Stack Overflow用户

发布于 2020-05-08 03:46:59

我不认为有一种完全矢量化的方法，只使用numpy函数(但我很高兴被证明是错误的！)

您可以通过使用itertools.accumulate来生成累积产品来隐藏循环。下面是一个例子。要创建A，我将使用使用scipy.stats.ortho_group生成的N随机正交矩阵来确保产品保持有界。

这里的前几行创建了形状(1000，4，4)的A。

In [101]: from scipy.stats import ortho_group

In [102]: N = 1000

In [103]: D = 4

In [104]: A = ortho_group.rvs(D, size=N)

使用itertools.accumulate计算累积产品，并将结果放入一个numpy数组中。

In [105]: from itertools import accumulate

In [106]: B = np.array(list(accumulate(A, np.matmul)))

验证我们从cumdot(A)获得了相同的结果。

In [107]: def cumdot(A): 
     ...:     B = np.empty(A.shape) 
     ...:     B[0] = A[0] 
     ...:     for i in range(1, A.shape[0]):   
     ...:         B[i] = B[i - 1] @ A[i] 
     ...:     return B 
     ...:                                                                                         

In [108]: B0 = cumdot(A)                                                                          

In [109]: (B == B0).all()                                                                         
Out[109]: True

检查一下表演。事实证明，使用itertools.accumulate的速度略快一些。

In [110]: %timeit B0 = cumdot(A)
2.89 ms ± 31.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [111]: %timeit B = np.array(list(accumulate(A, np.matmul)))
2.44 ms ± 33.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)