級聯矩陣乘法

Question

對於密集矩陣A和向量e，f的快速乘法A * diag（e）* A^T * f的任何建議。

這就是我現在擁有的。

v[:] = 0 
for i in range(N): 
    for j in range(N): 
     v[i] = v[i]+A[i,j]*e[j]*np.dot(A[:,j],f) 

感謝，

Answer 1

通過@rubenvb提出的建議可能是做到這一點的最簡單方法。另一種方法是使用einsum。

下面是一個例子。我將使用以下a，e和f：

In [95]: a 
Out[95]: 
array([[1, 2, 3], 
     [4, 5, 6], 
     [7, 8, 9]]) 

In [96]: e 
Out[96]: array([-1, 2, 3]) 

In [97]: f 
Out[97]: array([5, 4, 1]) 

這是直接翻譯你的公式爲numpy的代碼。這是基本相同@ rubenvb的建議：

In [98]: a.dot(np.diag(e)).dot(a.T).dot(f) 
Out[98]: array([ 556, 1132, 1708]) 

這裏的einsum版本：

In [99]: np.einsum('ij,j,jk,k', a, e, a.T, f) 
Out[99]: array([ 556, 1132, 1708]) 

可以不再需要通過交換與該參數關聯的索引標識移調a：

In [100]: np.einsum('ij,j,kj,k', a, e, a, f) 
Out[100]: array([ 556, 1132, 1708]) 

Answer 2

評論@ rubenvb's，它建議使用A.dot(np.diag(e)).dot(A.transpose()).dot(f)應該使它真的很快。但是，我們並不需要在diag(e)處設置2D數組，因此跳過了一個矩陣乘法。此外，我們可以交換A.T和f的位置，從而避免轉置。因此，簡化的和更有效的解決方案的發展，像這樣 -

A.dot(e*f.dot(A)) 

這裏的所有建議的方法體面大小的數組的快速運行測試 -

In [226]: # Setup inputs 
    ...: N = 200 
    ...: A = np.random.rand(N,N) 
    ...: e = np.random.rand(N,) 
    ...: f = np.random.rand(N,) 
    ...: 

In [227]: %timeit np.einsum('ij,j,kj,k', A, e, A, f) # @Warren Weckesser's soln 
10 loops, best of 3: 77.6 ms per loop 

In [228]: %timeit A.dot(np.diag(e)).dot(A.transpose()).dot(f) # @rubenvb's soln 
10 loops, best of 3: 18.6 ms per loop 

In [229]: %timeit A.dot(e*f.dot(A)) # Proposed here 
10000 loops, best of 3: 100 µs per loop