具有已知结构的矩阵的 NumPy 矩阵乘法效率

Question

我有两个要相乘的 NxN 矩阵：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(NNN) = 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 的低级实现细节可能会为此负责。所以我的问题是：如何利用矩阵 B 的结构来提高乘法效率而不牺牲 NumPy 的实现效率，在 c 中构建自己的低级矩阵乘法？

这种方法是对许多变量进行数值优化的一部分，因此，O(N^3) 是难以处理的，而 O(N^2) 可能会完成这项工作。

感谢您的帮助。另外，我是新手，所以请原谅任何新手错误。

Answer 1

我计时了，使用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 行数的增加，使用稀疏矩阵进行乘法的时间优势也应增加。

Answer 2

如果我正确理解了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