具有已知结构的矩阵的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(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次(对于必须乘以B的n列的A的每一行一次),必须发生N个标量乘法.

但是,第一个代码块(使用np.dot(A,B))要快得多.我知道(通过以下信息:为什么矩阵乘法比numpy更快,而不是Python中的ctypes？)np.dot的低级实现细节很可能归咎于此.所以我的问题是:如何在不牺牲NumPy的实现效率的情况下利用矩阵B的结构来提高乘法效率,而不在c中构建我自己的低级矩阵乘法？

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

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

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