python - 从两个稀疏矩阵计算任意行的点积的快速方法是什么

Question

例如...

import numpy as np
from scipy.sparse import csr_matrix

X = csr_matrix([[1,2,3], [4,5,6], [7,8,9]])
Y = csr_matrix([[1,2,3], [4,5,6], [7,8,9], [11,12,13]])

# Print matrices
X.toarray()
[[1, 2, 3],
 [4, 5, 6],
 [7, 8, 9]]

Y.toarray()
[[ 1,  2,  3],
 [ 4,  5,  6],
 [ 7,  8,  9],
 [11, 12, 13]]

我有一组索引对 (x,y)，分别代表 X 中的一行和 Y 中的一行。我想获取相应行的点积，但我不知道如何有效地做到这一点。

这是我尝试过的

# build arbitrary combinations of row from X and row from Y. Need to calculate dot product of each pair
x_idxs = np.array([2,2,1,0])
y_idxs = np.arange(Y.shape[0])

# current method (slow)
def get_dot_product(x_idx, y_idx):
    return np.dot(X[x_idx].toarray()[0], Y[y_idx].toarray()[0])

func_args = np.transpose(np.array([x_idxs, y_idxs]))
np.apply_along_axis(func1d=lambda x: get_dot_product(x[0], x[1]), axis=1, arr=func_args)

它可以工作，但随着X和Y变大而变慢。有没有更有效的方法？

更新

遵循 Warren 优雅但缓慢的解决方案，这里有一个更好的测试示例(以及基准)

X = csr_matrix(np.tile(np.repeat(1, 50000),(10000,1)))
Y = X
y_idxs = np.arange(Y.shape[0])
x_idxs = y_idxs

import time
start_time = time.time()
func_args = np.transpose(np.array([x_idxs, y_idxs]))
bg = np.apply_along_axis(func1d=lambda x: get_dot_product(x[0], x[1]), axis=1, arr=func_args)
print("--- %s seconds ---" % (time.time() - start_time)) # 15.48 seconds

start_time = time.time()
ww = X[x_idxs].multiply(Y[y_idxs]).sum(axis=1)
print("--- %s seconds ---" % (time.time() - start_time)) # 38.29 seconds

Answer 1

使用 X、Y、x_idxs 和 y_idxs，您可以执行以下操作:

In [160]: X[x_idxs].multiply(Y[y_idxs]).sum(axis=1)
Out[160]: 
matrix([[ 50],
        [122],
        [122],
        [ 74]])

它使用“奇特”索引(即使用任意序列进行索引以提取所需的行集)，然后进行逐点乘法和沿轴 1 求和来计算点积。

结果是一个 numpy 矩阵，您可以将其转换为常规 numpy 数组并根据需要进行展平。您甚至可以使用有点神秘的 A1 属性(getA1 方法的快捷方式):

In [178]: p = X[x_idxs].multiply(Y[y_idxs]).sum(axis=1)

In [179]: p
Out[179]: 
matrix([[ 50],
        [122],
        [122],
        [ 74]])

In [180]: p.A1
Out[180]: array([ 50, 122, 122,  74])

更新，有时间...

这是一个完整的脚本，用于将我的版本与原始版本的性能进行比较，使用实际上稀疏的数组X和Y(密度约为0.001，即大约0.1) %非零元素)。

import numpy as np
from scipy import sparse


def get_dot_product(x_idx, y_idx):
    return np.dot(X[x_idx].toarray()[0], Y[y_idx].toarray()[0])

print("Generating random sparse integer matrix X...")
X = (100000*sparse.rand(50000, 120000, density=0.001, format='csr')).astype(np.int64)
X.eliminate_zeros()
print("X has shape %s with %s nonzero elements." % (X.shape, X.nnz))
Y = X
y_idxs = np.arange(Y.shape[0])
x_idxs = y_idxs

import time
start_time = time.time()
func_args = np.transpose(np.array([x_idxs, y_idxs]))
bg = np.apply_along_axis(func1d=lambda x: get_dot_product(x[0], x[1]), axis=1, arr=func_args)
print("--- %8.5f seconds ---" % (time.time() - start_time))

start_time = time.time()
ww = X[x_idxs].multiply(Y[y_idxs]).sum(axis=1)
print("--- %8.5f seconds ---" % (time.time() - start_time))

输出:

Generating random sparse integer matrix X...
X has shape (50000, 120000) with 5999934 nonzero elements.
--- 18.29916 seconds ---
---  0.32749 seconds ---

对于不太稀疏的矩阵，速度差异不是很大，而对于足够稠密的矩阵，原始版本更快。