两个仅影响零值的稀疏矩阵的点积
Dot product of two sparse matrices affecting zero values only
我正在尝试计算一个简单的点积，但保留原始矩阵中的非零值不变。玩具示例：
import numpy as np
A = np.array([[2, 1, 1, 2],
[0, 2, 1, 0],
[1, 0, 1, 1],
[2, 2, 1, 0]])
B = np.array([[ 0.54331039, 0.41018682, 0.1582158 , 0.3486124 ],
[ 0.68804647, 0.29520239, 0.40654206, 0.20473451],
[ 0.69857579, 0.38958572, 0.30361365, 0.32256483],
[ 0.46195299, 0.79863505, 0.22431876, 0.59054473]])
期望的结果：
C = np.array([[ 2. , 1. , 1. , 2. ],
[ 2.07466874, 2. , 1. , 0.73203386],
[ 1. , 1.5984076 , 1. , 1. ],
[ 2. , 2. , 1. , 1.42925865]])
然而，所讨论的实际矩阵是稀疏的，看起来更像这样：
A = sparse.rand(250000, 1700, density=0.001, format='csr')
B = sparse.rand(1700, 1700, density=0.02, format='csr')
一个简单的方法就是使用掩码索引设置值，就像这样：
mask = A != 0
C = A.dot(B)
C[mask] = A[mask]
但是，我的原始数组很稀疏而且很大，因此通过索引赋值更改它们非常慢。转换为 lil 矩阵会有所帮助，但同样，转换本身也需要很多时间。
我想，另一种明显的方法就是求助于迭代并跳过屏蔽值，但我不想放弃 numpy/scipy-optimized 数组乘法的好处。
一些说明： 我实际上对某种特殊情况感兴趣，其中 B 总是正方形，因此，A 和 C形状相同。因此，如果有一个解决方案不适用于任意数组但适合我的情况，那很好。
更新： 一些尝试：
from scipy import sparse
import numpy as np
def naive(A, B):
mask = A != 0
out = A.dot(B).tolil()
out[mask] = A[mask]
return out.tocsr()
def proposed(A, B):
Az = A == 0
R, C = np.where(Az)
out = A.copy()
out[Az] = np.einsum('ij,ji->i', A[R], B[:, C])
return out
%timeit naive(A, B)
1 loops, best of 3: 4.04 s per loop
%timeit proposed(A, B)
/usr/local/lib/python2.7/dist-packages/scipy/sparse/compressed.py:215: SparseEfficiencyWarning: Comparing a sparse matrix with 0 using == is inefficient, try using != instead.
/usr/local/lib/python2.7/dist-packages/scipy/sparse/coo.pyc in __init__(self, arg1, shape, dtype, copy)
173 self.shape = M.shape
174
--> 175 self.row, self.col = M.nonzero()
176 self.data = M[self.row, self.col]
177 self.has_canonical_format = True
MemoryError:
另一项更新：
无法从 Cython 中获得更多或更少的有用信息，至少在与 Python 相去甚远的情况下是这样。想法是将点积留给 scipy 并尝试尽快设置这些原始值，如下所示：
cimport cython
@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef coo_replace(int [:] row1, int [:] col1, float [:] data1, int[:] row2, int[:] col2, float[:] data2):
cdef int N = row1.shape[0]
cdef int M = row2.shape[0]
cdef int i, j
cdef dict d = {}
for i in range(M):
d[(row2[i], col2[i])] = data2[i]
for j in range(N):
if (row1[j], col1[j]) in d:
data1[j] = d[(row1[j], col1[j])]
这比我之前的第一个 "naive" 实现（使用 .tolil()）好一点，但是按照 hpaulj 的方法，lil 可以被抛弃。也许用 std::map 之类的东西替换 python dict 会有所帮助。
破解了！好吧，我一路上学到了很多 scipy 特定于稀疏矩阵的东西。这是我可以召集的实现 -
# Find the indices in output array that are to be updated
R,C = ((A!=0).dot(B!=0)).nonzero()
mask = np.asarray(A[R,C]==0).ravel()
R,C = R[mask],C[mask]
# Make a copy of A and get the dot product through sliced rows and columns
# off A and B using the definition of matrix-multiplication
out = A.copy()
out[R,C] = (A[R].multiply(B[:,C].T).sum(1)).ravel()
最昂贵的部分似乎是逐元素乘法和求和。在一些快速计时测试中，这似乎对具有高度稀疏性的稀疏矩阵来说很好，可以在性能方面击败原始的基于点掩码的解决方案，我认为这来自于它对内存效率的关注。
运行时测试
函数定义-
def naive(A, B):
mask = A != 0
out = A.dot(B).tolil()
out[mask] = A[mask]
return out.tocsr()
def proposed(A, B):
R,C = ((A!=0).dot(B!=0)).nonzero()
mask = np.asarray(A[R,C]==0).ravel()
R,C = R[mask],C[mask]
out = A.copy()
out[R,C] = (A[R].multiply(B[:,C].T).sum(1)).ravel()
return out
计时 -
In [57]: # Input matrices
...: M,N = 25000, 170
...: A = sparse.rand(M, N, density=0.001, format='csr')
...: B = sparse.rand(N, N, density=0.02, format='csr')
...:
In [58]: %timeit naive(A, B)
10 loops, best of 3: 92.2 ms per loop
In [59]: %timeit proposed(A, B)
10 loops, best of 3: 132 ms per loop
In [60]: # Input matrices with increased sparse-ness
...: M,N = 25000, 170
...: A = sparse.rand(M, N, density=0.0001, format='csr')
...: B = sparse.rand(N, N, density=0.002, format='csr')
...:
In [61]: %timeit naive(A, B)
10 loops, best of 3: 78.1 ms per loop
In [62]: %timeit proposed(A, B)
100 loops, best of 3: 8.03 ms per loop
您的 naive 代码的一个可能更干净、更快的版本：
In [57]: r,c=A.nonzero() # this uses A.tocoo()
In [58]: C=A*B
In [59]: Cl=C.tolil()
In [60]: Cl[r,c]=A.tolil()[r,c]
In [61]: Cl.tocsr()
C[r,c]=A[r,c] 给出了效率警告，但我认为这更多的是为了让人们在循环中进行这种分配。
In [63]: %%timeit C=A*B
...: C[r,c]=A[r,c]
...
The slowest run took 7.32 times longer than the fastest....
1000 loops, best of 3: 334 µs per loop
In [64]: %%timeit C=A*B
...: Cl=C.tolil()
...: Cl[r,c]=A.tolil()[r,c]
...: Cl.tocsr()
...:
100 loops, best of 3: 2.83 ms per loop
我的 A 很小，只有 (250,100)，但看起来往返 lil 并不节省时间，尽管有警告。
当 A 稀疏
时，使用 A==0 进行掩码必然会出现问题
In [66]: Az=A==0
....SparseEfficiencyWarning...
In [67]: r1,c1=Az.nonzero()
与A的nonzeror相比，这个r1要大得多——稀疏矩阵中所有零的行索引；除了 25 个非零值以外的所有内容。
In [70]: r.shape
Out[70]: (25,)
In [71]: r1.shape
Out[71]: (24975,)
如果我用 r1 索引 A，我会得到一个更大的数组。实际上，我按其中的零数重复每一行
In [72]: A[r1,:]
Out[72]:
<24975x100 sparse matrix of type '<class 'numpy.float64'>'
with 2473 stored elements in Compressed Sparse Row format>
In [73]: A
Out[73]:
<250x100 sparse matrix of type '<class 'numpy.float64'>'
with 25 stored elements in Compressed Sparse Row format>
我将非零元素的形状和数量增加了大约 100（列数）。
定义 foo，并复制 Divakar 的测试：
def foo(A,B):
r,c = A.nonzero()
C = A*B
C[r,c] = A[r,c]
return C
In [83]: timeit naive(A,B)
100 loops, best of 3: 2.53 ms per loop
In [84]: timeit proposed(A,B)
/...
SparseEfficiencyWarning)
100 loops, best of 3: 4.48 ms per loop
In [85]: timeit foo(A,B)
...
SparseEfficiencyWarning)
100 loops, best of 3: 2.13 ms per loop
所以我的版本有适度的速度改进。正如 Divakar 发现的那样，改变稀疏度会改变相对优势。我希望尺寸也会改变它们。
A.nonzero 使用 coo 格式这一事实表明，使用该格式构建新数组可能是可行的。许多 sparse 代码通过 coo 值构建新矩阵。
In [97]: Co=C.tocoo()
In [98]: Ao=A.tocoo()
In [99]: r=np.concatenate((Co.row,Ao.row))
In [100]: c=np.concatenate((Co.col,Ao.col))
In [101]: d=np.concatenate((Co.data,Ao.data))
In [102]: r.shape
Out[102]: (79,)
In [103]: C1=sparse.csr_matrix((d,(r,c)),shape=A.shape)
In [104]: C1
Out[104]:
<250x100 sparse matrix of type '<class 'numpy.float64'>'
with 78 stored elements in Compressed Sparse Row format>
这个 C1 我认为与通过其他方式构造的 C 具有相同的非零元素。但我认为一个值是不同的，因为 r 更长。在这个特定的示例中，C 和 A 共享一个非零元素，并且 coo 样式的输入对这些元素求和，因为我们希望 A 值覆盖所有内容.
如果你能容忍这种差异，这是一种更快的方法（至少对于这个测试用例）：
def bar(A,B):
C=A*B
Co=C.tocoo()
Ao=A.tocoo()
r=np.concatenate((Co.row,Ao.row))
c=np.concatenate((Co.col,Ao.col))
d=np.concatenate((Co.data,Ao.data))
return sparse.csr_matrix((d,(r,c)),shape=A.shape)
In [107]: timeit bar(A,B)
1000 loops, best of 3: 1.03 ms per loop
Python 不是我的主要语言，但我认为这是一个有趣的问题，我想尝试一下 :)
预赛：
import numpy
import scipy.sparse
# example matrices and sparse versions
A = numpy.array([[1, 2, 0, 1], [1, 0, 1, 2], [0, 1, 2 ,1], [1, 2, 1, 0]])
B = numpy.array([[1,2,3,4],[1,2,3,4],[1,2,3,4],[1,2,3,4]])
A_s = scipy.sparse.lil_matrix(A)
B_s = scipy.sparse.lil_matrix(B)
原来你要转换的问题是：
C = A.dot(B)
C[A.nonzero()] = A[A.nonzero()]
稀疏的东西。
只是为了解决这个问题，上面的直接 "sparse" 翻译是：
C_s = A_s.dot(B_s)
C_s[A_s.nonzero()] = A_s[A_s.nonzero()]
但听起来你对此并不满意，因为它首先计算所有的点积，你担心这可能效率低下。
所以，你的问题是，如果你先找到零，然后只计算这些元素的点积，会更快吗？ IE。对于密集矩阵，这可能是这样的：
Xs, Ys = numpy.nonzero(A==0)
D = A[:]
D[Xs, Ys] = map ( lambda x,y: A[x,:].dot(B[:,y]), Xs, Ys)
让我们将其转换为稀疏矩阵。我在这里的主要绊脚石是找到 "Zero" 索引；由于 A_s==0 对稀疏矩阵没有意义，我是这样找到它们的：
Xmax, Ymax = A_s.shape
DenseSize = Xmax * Ymax
Xgrid, Ygrid = numpy.mgrid[0:Xmax, 0:Ymax]
Ygrid = Ygrid.reshape([DenseSize,1])[:,0]
Xgrid = Xgrid.reshape([DenseSize,1])[:,0]
AllIndices = numpy.array([Xgrid, Ygrid])
NonzeroIndices = numpy.array(A_s.nonzero())
ZeroIndices = numpy.array([x for x in AllIndices.T.tolist() if x not in NonzeroIndices.T.tolist()]).T
如果您知道更好/更快的方法，一定要尝试一下。一旦我们有了零索引，我们就可以像以前一样进行类似的映射：
D_s = A_s[:]
D_s[ZeroIndices[0], ZeroIndices[1]] = map ( lambda x, y : A_s[x,:].dot(B[:,y])[0], ZeroIndices[0], ZeroIndices[1] )
这会给出稀疏矩阵结果。
现在我不知道这是否更快。我试了一下主要是因为这是一个有趣的问题，看看我能否在 python 中完成。事实上，我怀疑它可能不会比直接全矩阵点积更快，因为它在大型数据集上使用列表推导和映射（就像你说的，你期望有很多零）。但它 是 对您 "how can I only calculate dot products for the zero values without doing multiplying the matrices as a whole" 问题的回答。我很想知道您是否尝试过它在数据集速度方面的比较。
编辑： 我在下面放了一个基于上述内容的示例 "block processing" 版本，我认为它应该允许您毫无问题地处理您的大型数据集。让我知道它是否有效。
import numpy
import scipy.sparse
# example matrices and sparse versions
A = numpy.array([[1, 2, 0, 1], [1, 0, 1, 2], [0, 1, 2 ,1], [1, 2, 1, 0]])
B = numpy.array([[1,2,3,4],[1,2,3,4],[1,2,3,4],[1,2,3,4]])
A_s = scipy.sparse.lil_matrix(A)
B_s = scipy.sparse.lil_matrix(B)
# Choose a grid division (i.e. how many processing blocks you want to create)
BlockGrid = numpy.array([2,2])
D_s = A_s[:] # initialise from A
Xmax, Ymax = A_s.shape
BaseBSiz = numpy.array([Xmax, Ymax]) / BlockGrid
for BIndX in range(0, Xmax, BlockGrid[0]):
for BIndY in range(0, Ymax, BlockGrid[1]):
BSizX, BSizY = D_s[ BIndX : BIndX + BaseBSiz[0], BIndY : BIndY + BaseBSiz[1] ].shape
Xgrid, Ygrid = numpy.mgrid[BIndX : BIndX + BSizX, BIndY : BIndY + BSizY]
Xgrid = Xgrid.reshape([BSizX*BSizY,1])[:,0]
Ygrid = Ygrid.reshape([BSizX*BSizY,1])[:,0]
AllInd = numpy.array([Xgrid, Ygrid]).T
NZeroInd = numpy.array(A_s[Xgrid, Ygrid].reshape((BSizX,BSizY)).nonzero()).T + numpy.array([[BIndX],[BIndY]]).T
ZeroInd = numpy.array([x for x in AllInd.tolist() if x not in NZeroInd.tolist()]).T
#
# Replace zero-values in current block
D_s[ZeroInd[0], ZeroInd[1]] = map ( lambda x, y : A_s[x,:].dot(B[:,y])[0], ZeroInd[0], ZeroInd[1] )
我正在尝试计算一个简单的点积，但保留原始矩阵中的非零值不变。玩具示例：
import numpy as np
A = np.array([[2, 1, 1, 2],
[0, 2, 1, 0],
[1, 0, 1, 1],
[2, 2, 1, 0]])
B = np.array([[ 0.54331039, 0.41018682, 0.1582158 , 0.3486124 ],
[ 0.68804647, 0.29520239, 0.40654206, 0.20473451],
[ 0.69857579, 0.38958572, 0.30361365, 0.32256483],
[ 0.46195299, 0.79863505, 0.22431876, 0.59054473]])
期望的结果：
C = np.array([[ 2. , 1. , 1. , 2. ],
[ 2.07466874, 2. , 1. , 0.73203386],
[ 1. , 1.5984076 , 1. , 1. ],
[ 2. , 2. , 1. , 1.42925865]])
然而，所讨论的实际矩阵是稀疏的，看起来更像这样：
A = sparse.rand(250000, 1700, density=0.001, format='csr')
B = sparse.rand(1700, 1700, density=0.02, format='csr')
一个简单的方法就是使用掩码索引设置值，就像这样：
mask = A != 0
C = A.dot(B)
C[mask] = A[mask]
但是，我的原始数组很稀疏而且很大，因此通过索引赋值更改它们非常慢。转换为 lil 矩阵会有所帮助，但同样，转换本身也需要很多时间。
我想，另一种明显的方法就是求助于迭代并跳过屏蔽值，但我不想放弃 numpy/scipy-optimized 数组乘法的好处。
一些说明： 我实际上对某种特殊情况感兴趣，其中 B 总是正方形，因此，A 和 C形状相同。因此，如果有一个解决方案不适用于任意数组但适合我的情况，那很好。
更新： 一些尝试：
from scipy import sparse
import numpy as np
def naive(A, B):
mask = A != 0
out = A.dot(B).tolil()
out[mask] = A[mask]
return out.tocsr()
def proposed(A, B):
Az = A == 0
R, C = np.where(Az)
out = A.copy()
out[Az] = np.einsum('ij,ji->i', A[R], B[:, C])
return out
%timeit naive(A, B)
1 loops, best of 3: 4.04 s per loop
%timeit proposed(A, B)
/usr/local/lib/python2.7/dist-packages/scipy/sparse/compressed.py:215: SparseEfficiencyWarning: Comparing a sparse matrix with 0 using == is inefficient, try using != instead.
/usr/local/lib/python2.7/dist-packages/scipy/sparse/coo.pyc in __init__(self, arg1, shape, dtype, copy)
173 self.shape = M.shape
174
--> 175 self.row, self.col = M.nonzero()
176 self.data = M[self.row, self.col]
177 self.has_canonical_format = True
MemoryError:
另一项更新：
无法从 Cython 中获得更多或更少的有用信息，至少在与 Python 相去甚远的情况下是这样。想法是将点积留给 scipy 并尝试尽快设置这些原始值，如下所示：
cimport cython
@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef coo_replace(int [:] row1, int [:] col1, float [:] data1, int[:] row2, int[:] col2, float[:] data2):
cdef int N = row1.shape[0]
cdef int M = row2.shape[0]
cdef int i, j
cdef dict d = {}
for i in range(M):
d[(row2[i], col2[i])] = data2[i]
for j in range(N):
if (row1[j], col1[j]) in d:
data1[j] = d[(row1[j], col1[j])]
这比我之前的第一个 "naive" 实现（使用 .tolil()）好一点，但是按照 hpaulj 的方法，lil 可以被抛弃。也许用 std::map 之类的东西替换 python dict 会有所帮助。
破解了！好吧，我一路上学到了很多 scipy 特定于稀疏矩阵的东西。这是我可以召集的实现 -
# Find the indices in output array that are to be updated
R,C = ((A!=0).dot(B!=0)).nonzero()
mask = np.asarray(A[R,C]==0).ravel()
R,C = R[mask],C[mask]
# Make a copy of A and get the dot product through sliced rows and columns
# off A and B using the definition of matrix-multiplication
out = A.copy()
out[R,C] = (A[R].multiply(B[:,C].T).sum(1)).ravel()
最昂贵的部分似乎是逐元素乘法和求和。在一些快速计时测试中，这似乎对具有高度稀疏性的稀疏矩阵来说很好，可以在性能方面击败原始的基于点掩码的解决方案，我认为这来自于它对内存效率的关注。
运行时测试
函数定义-
def naive(A, B):
mask = A != 0
out = A.dot(B).tolil()
out[mask] = A[mask]
return out.tocsr()
def proposed(A, B):
R,C = ((A!=0).dot(B!=0)).nonzero()
mask = np.asarray(A[R,C]==0).ravel()
R,C = R[mask],C[mask]
out = A.copy()
out[R,C] = (A[R].multiply(B[:,C].T).sum(1)).ravel()
return out
计时 -
In [57]: # Input matrices
...: M,N = 25000, 170
...: A = sparse.rand(M, N, density=0.001, format='csr')
...: B = sparse.rand(N, N, density=0.02, format='csr')
...:
In [58]: %timeit naive(A, B)
10 loops, best of 3: 92.2 ms per loop
In [59]: %timeit proposed(A, B)
10 loops, best of 3: 132 ms per loop
In [60]: # Input matrices with increased sparse-ness
...: M,N = 25000, 170
...: A = sparse.rand(M, N, density=0.0001, format='csr')
...: B = sparse.rand(N, N, density=0.002, format='csr')
...:
In [61]: %timeit naive(A, B)
10 loops, best of 3: 78.1 ms per loop
In [62]: %timeit proposed(A, B)
100 loops, best of 3: 8.03 ms per loop
您的 naive 代码的一个可能更干净、更快的版本：
In [57]: r,c=A.nonzero() # this uses A.tocoo()
In [58]: C=A*B
In [59]: Cl=C.tolil()
In [60]: Cl[r,c]=A.tolil()[r,c]
In [61]: Cl.tocsr()
C[r,c]=A[r,c] 给出了效率警告，但我认为这更多的是为了让人们在循环中进行这种分配。
In [63]: %%timeit C=A*B
...: C[r,c]=A[r,c]
...
The slowest run took 7.32 times longer than the fastest....
1000 loops, best of 3: 334 µs per loop
In [64]: %%timeit C=A*B
...: Cl=C.tolil()
...: Cl[r,c]=A.tolil()[r,c]
...: Cl.tocsr()
...:
100 loops, best of 3: 2.83 ms per loop
我的 A 很小，只有 (250,100)，但看起来往返 lil 并不节省时间，尽管有警告。
当 A 稀疏
A==0 进行掩码必然会出现问题
In [66]: Az=A==0
....SparseEfficiencyWarning...
In [67]: r1,c1=Az.nonzero()
与A的nonzeror相比，这个r1要大得多——稀疏矩阵中所有零的行索引；除了 25 个非零值以外的所有内容。
In [70]: r.shape
Out[70]: (25,)
In [71]: r1.shape
Out[71]: (24975,)
如果我用 r1 索引 A，我会得到一个更大的数组。实际上，我按其中的零数重复每一行
In [72]: A[r1,:]
Out[72]:
<24975x100 sparse matrix of type '<class 'numpy.float64'>'
with 2473 stored elements in Compressed Sparse Row format>
In [73]: A
Out[73]:
<250x100 sparse matrix of type '<class 'numpy.float64'>'
with 25 stored elements in Compressed Sparse Row format>
我将非零元素的形状和数量增加了大约 100（列数）。
定义 foo，并复制 Divakar 的测试：
def foo(A,B):
r,c = A.nonzero()
C = A*B
C[r,c] = A[r,c]
return C
In [83]: timeit naive(A,B)
100 loops, best of 3: 2.53 ms per loop
In [84]: timeit proposed(A,B)
/...
SparseEfficiencyWarning)
100 loops, best of 3: 4.48 ms per loop
In [85]: timeit foo(A,B)
...
SparseEfficiencyWarning)
100 loops, best of 3: 2.13 ms per loop
所以我的版本有适度的速度改进。正如 Divakar 发现的那样，改变稀疏度会改变相对优势。我希望尺寸也会改变它们。
A.nonzero 使用 coo 格式这一事实表明，使用该格式构建新数组可能是可行的。许多 sparse 代码通过 coo 值构建新矩阵。
In [97]: Co=C.tocoo()
In [98]: Ao=A.tocoo()
In [99]: r=np.concatenate((Co.row,Ao.row))
In [100]: c=np.concatenate((Co.col,Ao.col))
In [101]: d=np.concatenate((Co.data,Ao.data))
In [102]: r.shape
Out[102]: (79,)
In [103]: C1=sparse.csr_matrix((d,(r,c)),shape=A.shape)
In [104]: C1
Out[104]:
<250x100 sparse matrix of type '<class 'numpy.float64'>'
with 78 stored elements in Compressed Sparse Row format>
这个 C1 我认为与通过其他方式构造的 C 具有相同的非零元素。但我认为一个值是不同的，因为 r 更长。在这个特定的示例中，C 和 A 共享一个非零元素，并且 coo 样式的输入对这些元素求和，因为我们希望 A 值覆盖所有内容.
如果你能容忍这种差异，这是一种更快的方法（至少对于这个测试用例）：
def bar(A,B):
C=A*B
Co=C.tocoo()
Ao=A.tocoo()
r=np.concatenate((Co.row,Ao.row))
c=np.concatenate((Co.col,Ao.col))
d=np.concatenate((Co.data,Ao.data))
return sparse.csr_matrix((d,(r,c)),shape=A.shape)
In [107]: timeit bar(A,B)
1000 loops, best of 3: 1.03 ms per loop
Python 不是我的主要语言，但我认为这是一个有趣的问题，我想尝试一下 :)
预赛：
import numpy
import scipy.sparse
# example matrices and sparse versions
A = numpy.array([[1, 2, 0, 1], [1, 0, 1, 2], [0, 1, 2 ,1], [1, 2, 1, 0]])
B = numpy.array([[1,2,3,4],[1,2,3,4],[1,2,3,4],[1,2,3,4]])
A_s = scipy.sparse.lil_matrix(A)
B_s = scipy.sparse.lil_matrix(B)
原来你要转换的问题是：
C = A.dot(B)
C[A.nonzero()] = A[A.nonzero()]
稀疏的东西。 只是为了解决这个问题，上面的直接 "sparse" 翻译是：
C_s = A_s.dot(B_s)
C_s[A_s.nonzero()] = A_s[A_s.nonzero()]
但听起来你对此并不满意，因为它首先计算所有的点积，你担心这可能效率低下。
所以，你的问题是，如果你先找到零，然后只计算这些元素的点积，会更快吗？ IE。对于密集矩阵，这可能是这样的：
Xs, Ys = numpy.nonzero(A==0)
D = A[:]
D[Xs, Ys] = map ( lambda x,y: A[x,:].dot(B[:,y]), Xs, Ys)
让我们将其转换为稀疏矩阵。我在这里的主要绊脚石是找到 "Zero" 索引；由于 A_s==0 对稀疏矩阵没有意义，我是这样找到它们的：
Xmax, Ymax = A_s.shape
DenseSize = Xmax * Ymax
Xgrid, Ygrid = numpy.mgrid[0:Xmax, 0:Ymax]
Ygrid = Ygrid.reshape([DenseSize,1])[:,0]
Xgrid = Xgrid.reshape([DenseSize,1])[:,0]
AllIndices = numpy.array([Xgrid, Ygrid])
NonzeroIndices = numpy.array(A_s.nonzero())
ZeroIndices = numpy.array([x for x in AllIndices.T.tolist() if x not in NonzeroIndices.T.tolist()]).T
如果您知道更好/更快的方法，一定要尝试一下。一旦我们有了零索引，我们就可以像以前一样进行类似的映射：
D_s = A_s[:]
D_s[ZeroIndices[0], ZeroIndices[1]] = map ( lambda x, y : A_s[x,:].dot(B[:,y])[0], ZeroIndices[0], ZeroIndices[1] )
这会给出稀疏矩阵结果。
现在我不知道这是否更快。我试了一下主要是因为这是一个有趣的问题，看看我能否在 python 中完成。事实上，我怀疑它可能不会比直接全矩阵点积更快，因为它在大型数据集上使用列表推导和映射（就像你说的，你期望有很多零）。但它 是 对您 "how can I only calculate dot products for the zero values without doing multiplying the matrices as a whole" 问题的回答。我很想知道您是否尝试过它在数据集速度方面的比较。
编辑： 我在下面放了一个基于上述内容的示例 "block processing" 版本，我认为它应该允许您毫无问题地处理您的大型数据集。让我知道它是否有效。
import numpy
import scipy.sparse
# example matrices and sparse versions
A = numpy.array([[1, 2, 0, 1], [1, 0, 1, 2], [0, 1, 2 ,1], [1, 2, 1, 0]])
B = numpy.array([[1,2,3,4],[1,2,3,4],[1,2,3,4],[1,2,3,4]])
A_s = scipy.sparse.lil_matrix(A)
B_s = scipy.sparse.lil_matrix(B)
# Choose a grid division (i.e. how many processing blocks you want to create)
BlockGrid = numpy.array([2,2])
D_s = A_s[:] # initialise from A
Xmax, Ymax = A_s.shape
BaseBSiz = numpy.array([Xmax, Ymax]) / BlockGrid
for BIndX in range(0, Xmax, BlockGrid[0]):
for BIndY in range(0, Ymax, BlockGrid[1]):
BSizX, BSizY = D_s[ BIndX : BIndX + BaseBSiz[0], BIndY : BIndY + BaseBSiz[1] ].shape
Xgrid, Ygrid = numpy.mgrid[BIndX : BIndX + BSizX, BIndY : BIndY + BSizY]
Xgrid = Xgrid.reshape([BSizX*BSizY,1])[:,0]
Ygrid = Ygrid.reshape([BSizX*BSizY,1])[:,0]
AllInd = numpy.array([Xgrid, Ygrid]).T
NZeroInd = numpy.array(A_s[Xgrid, Ygrid].reshape((BSizX,BSizY)).nonzero()).T + numpy.array([[BIndX],[BIndY]]).T
ZeroInd = numpy.array([x for x in AllInd.tolist() if x not in NZeroInd.tolist()]).T
#
# Replace zero-values in current block
D_s[ZeroInd[0], ZeroInd[1]] = map ( lambda x, y : A_s[x,:].dot(B[:,y])[0], ZeroInd[0], ZeroInd[1] )