NumPy的矩陣乘法效率的矩陣結構已知

Question

我有我要乘在一起的兩個N×N的矩陣：A和B.在NumPy的，我用：

import numpy as np 
C = np.dot(A, B) 

但是，我碰巧知道矩陣B只有第n行和第n列非零（這直接來自產生矩陣的分析公式，毫無疑問總是這種情況）。

希望能利用這一事實的優點和減少以生成C需要的乘法的數量，我取代了以上：

import numpy as np 
for row in range(0, N): 
    for col in range(0, N): 
     if col != n: 
      C[row, col] = A[row, n]*B[n, col] #Just one scalar multiplication 
     else: 
      C[row, col] = np.dot(A[row], B[:, n]) 

分析上，這應該減少總複雜性如下：在一般的情況下（不使用任何花哨的技巧，只是基本的矩陣乘法）C = AB，其中A和B都是NxN，應該是O（N^3）。也就是說，所有N行必須乘以所有N列，並且這些點積中的每一個都包含N次乘法=> O（N N N）= O（N^3）＃

利用B的結構作爲我做過上面的事情，但是應該是O（N^2 + N^2）= O（2N^2）= O（N^2）。也就是說，所有N行必須乘以所有N列，但是，對於所有這些（除了涉及'B [：，n]'的那些行），只需要一個標量乘法：只有'B [:, m]'中的一個元素對於m！= n是非零的。當n == m時，會發生N次（每次A行必須乘以B的列n），必須發生N次標量乘法。＃

但是，第一塊代碼（使用np.dot （A，B））實質上更快。我知道（通過如下信息：Why is matrix multiplication faster with numpy than with ctypes in Python?）np.dot的低層實現細節很可能會歸咎於此。所以我的問題是：如何在不犧牲NumPy的實現效率的情況下利用矩陣B的結構來提高乘法效率，而無需在c中構建自己的低級矩陣乘法？

該方法是對許多變量進行數值優化的一部分，因此，O（N^3）是棘手的，而O（N^2）可能會完成工作。

謝謝你的幫助。另外，我是新手，所以請原諒任何新手錯誤。

Answer 1

如果我的理解A和B正確的，那麼我不明白for循環，爲什麼你不只是由兩個非零向量相乘：

# say A & B are like this: 
n, N = 3, 5 
A = np.array(np.random.randn(N, N)) 

B = np.zeros_like(A) 
B[ n ] = np.random.randn(N) 
B[:, n] = np.random.randn(N) 

採取非零行和列B的：

rowb, colb = B[n,:], np.copy(B[:,n]) 
colb[ n ] = 0 

乘法A由這兩個向量：

X = np.outer(A[:,n], rowb) 
X[:,n] += np.dot(A, colb) 

驗證檢查：

X - np.dot(A, B) 

與N=100：

%timeit np.dot(A, B) 
1000 loops, best of 3: 1.39 ms per loop 

%timeit colb = np.copy(B[:,n]); colb[ n ] = 0; X = np.outer(A[:,n], B[n,:]); X[:,n] += np.dot(A, colb) 
10000 loops, best of 3: 98.5 µs per loop 

Answer 2

我計時，以及如何使用sparse更快：

import numpy as np 
from scipy import sparse 

from timeit import timeit 

A = np.random.rand(100,100) 
B = np.zeros(A.shape, dtype=np.float) 

B[3] = np.random.rand(100) 
B[:,3] = np.random.rand(100) 

sparse_B = sparse.csr_matrix(B) 

n = 1000 

t1 = timeit('np.dot(A, B)', 'from __main__ import np, A, B', number=n) 
print 'dense way : {}'.format(t1) 
t2 = timeit('A * sparse_B', 'from __main__ import A, sparse_B',number=n) 
print 'sparse way : {}'.format(t2) 

結果：

dense way : 1.15117192268 
sparse way : 0.113152980804 
>>> np.allclose(np.dot(A, B), A * sparse_B) 
True 

隨着B的行數增加，使用稀疏矩陣乘法的時間優勢也應該增加。