Numpy/Scipy 广播计算某个元素的标量积
Numpy/Scipy broadcast calculating scalar product for a certain elements
我有一个像 A 这样的稀疏矩阵
和一个数据框 (df)，其中包含用于计算标量积的行。
Row1 Row2 Value
2 147 scalar product of vectors at Row1 and Raw2 in matrix A
我可以在不循环等情况下以广播方式进行吗？
在我的例子中，A 大小为 1m*100k，数据帧为 10M
如果我没有正确理解你的问题，你可以使用 Pandas 中的 dot 函数来计算两个系列之间的点积：
A['Row1'].dot(A['Row2'])
文档：
http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.dot.html
我猜.assign()和.apply()（对于pandas > 0.16.0）是合适的：
import numpy as np
from pandas import DataFrame
from scipy.sparse import bsr_matrix
df = DataFrame(np.random.randint(4, size=(4, 2)), columns=['Row1', 'Row2'])
A = bsr_matrix([[1, 2, 3],
[2, 1, 4],
[0, 2, 2],
[3, 0, 3]])
A = A.tocsr() # Skip this if your matrix is csc_, csr_, dok_ or lil_matrix
df.assign(Value=df.apply(lambda row: A[row[0]].dot(A[row[1]].transpose())[0, 0], axis=1))
Out[15]:
Row1 Row2 Value
0 1 3 18
1 1 0 16
2 0 0 14
3 3 2 6
从一个小 'sparse' 矩阵开始（csr 最适合数学）：
In [167]: A=sparse.csr_matrix([[1, 2, 3], # Vadim's example
[2, 1, 4],
[0, 2, 2],
[3, 0, 3]])
In [168]: AA=A.A # dense equivalent
In [169]: idx=np.array([[1,1,0,3],[3,0,0,2]]).T # indexes
我会坚持使用 numpy 版本（Pandas 构建在 numpy 之上）
我们可以采用所有行点积，select 由 idx:
定义的子集
In [170]: (AA.dot(AA.T))[idx[:,0], idx[:,1]]
Out[170]: array([18, 16, 14, 6], dtype=int32)
稀疏矩阵乘积（A.dot(A.T) 也适用：
In [171]: (A*A.T)[idx[:,0], idx[:,1]]
Out[171]: matrix([[18, 16, 14, 6]], dtype=int32)
或者我们可以先select行，然后取乘积之和。我们不想在这里使用 dot，因为我们没有采用所有组合。
In [172]: (AA[idx[:,0]]*AA[idx[:,1]]).sum(axis=1)
Out[172]: array([18, 16, 14, 6], dtype=int32)
此计算的 einsum 版本：
In [180]: np.einsum('ij,ij->i',AA[idx[:,0]],AA[idx[:,1]])
Out[180]: array([18, 16, 14, 6], dtype=int32)
sparse也可以这样做（*是矩阵乘积，.multiply是逐个元素）。
In [173]: (A[idx[:,0]].multiply(A[idx[:,1]])).sum(axis=1)
Out[173]:
matrix([[18],
[16],
[14],
[ 6]], dtype=int32)
对于这个小案例，密集版本更快。稀疏行索引很慢。
In [181]: timeit np.einsum('ij,ij->i', AA[idx[:,0]], AA[idx[:,1]])
100000 loops, best of 3: 18.1 µs per loop
In [182]: timeit (A[idx[:,0]].multiply(A[idx[:,1]])).sum(axis=1)
1000 loops, best of 3: 1.32 ms per loop
In [184]: timeit (AA.dot(AA.T))[idx[:,0], idx[:,1]]
100000 loops, best of 3: 9.62 µs per loop
In [185]: timeit (A*A.T)[idx[:,0], idx[:,1]]
1000 loops, best of 3: 689 µs per loop
我差点忘了 - 迭代版本：
In [191]: timeit [AA[i].dot(AA[j]) for i,j in idx]
10000 loops, best of 3: 38.4 µs per loop
In [192]: timeit [A[i].multiply(A[j]).sum() for i,j in idx]
100 loops, best of 3: 2.58 ms per loop
lil 格式矩阵的行索引更快
In [207]: Al=A.tolil()
In [208]: timeit A[idx[:,0]]
1000 loops, best of 3: 476 µs per loop
In [209]: timeit Al[idx[:,0]]
1000 loops, best of 3: 234 µs per loop
但是当它转换回 csr 进行乘法运算时，它可能不会节省时间。
===============
在其他最近的 SO 问题中，我讨论了更快地索引稀疏行或列的方法。但在这些中，最终目标是对一组 select 行或列求和。为此，使用矩阵乘积实际上是最快的——使用 1 和 0 的矩阵。在这里应用这个想法有点棘手。
查看 csr.__getitem__ 索引函数，我发现它实际上使用矩阵乘积进行 A[idx,:] 索引。它创建一个 extractor 矩阵，其功能如下：
def extractor(indices,N):
"""Return a sparse matrix P so that P*self implements
slicing of the form self[[1,2,3],:]
"""
indptr = np.arange(len(indices)+1, dtype=int)
data = np.ones(len(indices), dtype=int)
shape = (len(indices),N)
return sparse.csr_matrix((data,indices,indptr), shape=shape)
In [328]: %%timeit
.....: A1=extractor(idx[:,0],4)*A
.....: A2=extractor(idx[:,1],4)*A
.....: (A1.multiply(A2)).sum(axis=1)
.....:
1000 loops, best of 3: 1.14 ms per loop
这次比用 A[idx[:,0],:]（上面的 In[182]）产生的稍好 - 大概是因为它稍微简化了操作。它应该以相同的方式缩放。
之所以有效，是因为 idx0 是从 [1,1,0,3]
派生的布尔矩阵
In [330]: extractor(idx[:,0],4).A
Out[330]:
array([[0, 1, 0, 0],
[0, 1, 0, 0],
[1, 0, 0, 0],
[0, 0, 0, 1]])
In [296]: A[idx[:,0],:].A
Out[296]:
array([[2, 1, 4],
[2, 1, 4],
[1, 2, 3],
[3, 0, 3]], dtype=int32)
In [331]: (extractor(idx[:,0],4)*A).A
Out[331]:
array([[2, 1, 4],
[2, 1, 4],
[1, 2, 3],
[3, 0, 3]], dtype=int32)
================
总而言之，如果问题太大而无法直接使用密集数组，那么扩展到大型稀疏情况的最佳选择是
(A[idx[:,0]].multiply(A[idx[:,1]])).sum(axis=1)
如果这仍然产生内存错误，则迭代，可能遍历 idx 数组（或数据帧）的组。
我有一个像 A 这样的稀疏矩阵
和一个数据框 (df)，其中包含用于计算标量积的行。
Row1 Row2 Value
2 147 scalar product of vectors at Row1 and Raw2 in matrix A
我可以在不循环等情况下以广播方式进行吗？
在我的例子中，A 大小为 1m*100k，数据帧为 10M
如果我没有正确理解你的问题，你可以使用 Pandas 中的 dot 函数来计算两个系列之间的点积：
A['Row1'].dot(A['Row2'])
文档： http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.dot.html
我猜.assign()和.apply()（对于pandas > 0.16.0）是合适的：
import numpy as np
from pandas import DataFrame
from scipy.sparse import bsr_matrix
df = DataFrame(np.random.randint(4, size=(4, 2)), columns=['Row1', 'Row2'])
A = bsr_matrix([[1, 2, 3],
[2, 1, 4],
[0, 2, 2],
[3, 0, 3]])
A = A.tocsr() # Skip this if your matrix is csc_, csr_, dok_ or lil_matrix
df.assign(Value=df.apply(lambda row: A[row[0]].dot(A[row[1]].transpose())[0, 0], axis=1))
Out[15]:
Row1 Row2 Value
0 1 3 18
1 1 0 16
2 0 0 14
3 3 2 6
从一个小 'sparse' 矩阵开始（csr 最适合数学）：
In [167]: A=sparse.csr_matrix([[1, 2, 3], # Vadim's example
[2, 1, 4],
[0, 2, 2],
[3, 0, 3]])
In [168]: AA=A.A # dense equivalent
In [169]: idx=np.array([[1,1,0,3],[3,0,0,2]]).T # indexes
我会坚持使用 numpy 版本（Pandas 构建在 numpy 之上）
我们可以采用所有行点积，select 由 idx:
In [170]: (AA.dot(AA.T))[idx[:,0], idx[:,1]]
Out[170]: array([18, 16, 14, 6], dtype=int32)
稀疏矩阵乘积（A.dot(A.T) 也适用：
In [171]: (A*A.T)[idx[:,0], idx[:,1]]
Out[171]: matrix([[18, 16, 14, 6]], dtype=int32)
或者我们可以先select行，然后取乘积之和。我们不想在这里使用 dot，因为我们没有采用所有组合。
In [172]: (AA[idx[:,0]]*AA[idx[:,1]]).sum(axis=1)
Out[172]: array([18, 16, 14, 6], dtype=int32)
此计算的 einsum 版本：
In [180]: np.einsum('ij,ij->i',AA[idx[:,0]],AA[idx[:,1]])
Out[180]: array([18, 16, 14, 6], dtype=int32)
sparse也可以这样做（*是矩阵乘积，.multiply是逐个元素）。
In [173]: (A[idx[:,0]].multiply(A[idx[:,1]])).sum(axis=1)
Out[173]:
matrix([[18],
[16],
[14],
[ 6]], dtype=int32)
对于这个小案例，密集版本更快。稀疏行索引很慢。
In [181]: timeit np.einsum('ij,ij->i', AA[idx[:,0]], AA[idx[:,1]])
100000 loops, best of 3: 18.1 µs per loop
In [182]: timeit (A[idx[:,0]].multiply(A[idx[:,1]])).sum(axis=1)
1000 loops, best of 3: 1.32 ms per loop
In [184]: timeit (AA.dot(AA.T))[idx[:,0], idx[:,1]]
100000 loops, best of 3: 9.62 µs per loop
In [185]: timeit (A*A.T)[idx[:,0], idx[:,1]]
1000 loops, best of 3: 689 µs per loop
我差点忘了 - 迭代版本：
In [191]: timeit [AA[i].dot(AA[j]) for i,j in idx]
10000 loops, best of 3: 38.4 µs per loop
In [192]: timeit [A[i].multiply(A[j]).sum() for i,j in idx]
100 loops, best of 3: 2.58 ms per loop
lil 格式矩阵的行索引更快
In [207]: Al=A.tolil()
In [208]: timeit A[idx[:,0]]
1000 loops, best of 3: 476 µs per loop
In [209]: timeit Al[idx[:,0]]
1000 loops, best of 3: 234 µs per loop
但是当它转换回 csr 进行乘法运算时，它可能不会节省时间。
===============
在其他最近的 SO 问题中，我讨论了更快地索引稀疏行或列的方法。但在这些中，最终目标是对一组 select 行或列求和。为此，使用矩阵乘积实际上是最快的——使用 1 和 0 的矩阵。在这里应用这个想法有点棘手。
查看 csr.__getitem__ 索引函数，我发现它实际上使用矩阵乘积进行 A[idx,:] 索引。它创建一个 extractor 矩阵，其功能如下：
def extractor(indices,N):
"""Return a sparse matrix P so that P*self implements
slicing of the form self[[1,2,3],:]
"""
indptr = np.arange(len(indices)+1, dtype=int)
data = np.ones(len(indices), dtype=int)
shape = (len(indices),N)
return sparse.csr_matrix((data,indices,indptr), shape=shape)
In [328]: %%timeit
.....: A1=extractor(idx[:,0],4)*A
.....: A2=extractor(idx[:,1],4)*A
.....: (A1.multiply(A2)).sum(axis=1)
.....:
1000 loops, best of 3: 1.14 ms per loop
这次比用 A[idx[:,0],:]（上面的 In[182]）产生的稍好 - 大概是因为它稍微简化了操作。它应该以相同的方式缩放。
之所以有效，是因为 idx0 是从 [1,1,0,3]
In [330]: extractor(idx[:,0],4).A
Out[330]:
array([[0, 1, 0, 0],
[0, 1, 0, 0],
[1, 0, 0, 0],
[0, 0, 0, 1]])
In [296]: A[idx[:,0],:].A
Out[296]:
array([[2, 1, 4],
[2, 1, 4],
[1, 2, 3],
[3, 0, 3]], dtype=int32)
In [331]: (extractor(idx[:,0],4)*A).A
Out[331]:
array([[2, 1, 4],
[2, 1, 4],
[1, 2, 3],
[3, 0, 3]], dtype=int32)
================
总而言之，如果问题太大而无法直接使用密集数组，那么扩展到大型稀疏情况的最佳选择是
(A[idx[:,0]].multiply(A[idx[:,1]])).sum(axis=1)
如果这仍然产生内存错误，则迭代，可能遍历 idx 数组（或数据帧）的组。