假设我有一个稀疏数组和一个密集数组，它们的列数相同但行数较少:
from scipy.sparse import csr_matrix
import numpy as np
sp_arr = csr_matrix(np.array([[1,0,0,0,1],[0,0,1,0,0],[0,1,0,0,1],[0,0,0,1,1],[0,0,0,1,0]]))
arr = np.random.rand(10).reshape(2,5)
print(arr)
[[ 0.47027789 0.82510323 0.01321617 0.66264852 0.3618022 ]
[ 0.80198907 0.36350616 0.10254934 0.65209401 0.094961 ]]
我想得到一个数组，其中包含索引的所有子矩阵，这些子矩阵包含稀疏数组每一行的值。例如sp_arr中数据的索引如下:
0: [0, 4] 1: [2] 2: [1, 4] 3: [3, 4] 4: [3]
我的输出应该是这样的:
array([array([[ 0.47027789, 0.3618022 ],
[ 0.80198907, 0.094961 ]]),
array([[ 0.01321617],
[ 0.10254934]]),
array([[ 0.82510323, 0.3618022 ],
[ 0.36350616, 0.094961 ]]),
array([[ 0.66264852, 0.3618022 ],
[ 0.65209401, 0.094961 ]]),
array([[ 0.66264852],
[ 0.65209401]])], dtype=object)
我可以使用以下代码创建它，但是随着数组大小的增加(在我的情况下很大)，它变得非常慢。
output = np.empty(sp_arr.shape[0], dtype=object)
for row in range(sp_arr.shape[0]):
output[row] = arr[:, sp_arr[row].indices]
我考虑过将过程矢量化并沿轴应用它，但是 np.apply_along_axis 不适用于稀疏矩阵，不幸的是，虽然这个例子足够小，可以使密集然后使用 apply_along_axis 我的实际稀疏矩阵太大而无法这样做 (>100Gb)。
我曾想也许有一种奇特的方法来索引或使用 hsplit 之类的东西来使用已经矢量化的方法来实现这一点，但到目前为止我还没有运气。有什么办法可以让我逃避吗？
更新
根据@Divakar 的回答，这很棒，我找到了另一种方法来实现同样的事情，并且改进非常微小，可以忽略不计。
@Divakars 最佳答案是:
def app2(sp_arr, arr):
r,c = sp_arr.nonzero()
idx = np.flatnonzero(r[1:] > r[:-1])+1
idx0 = np.concatenate(( [0] , idx, [r.size] ))
arr_c = arr[:,c]
return [arr_c[:,i:j] for i,j in zip(idx0[:-1], idx0[1:])]
这让我的表现提高了 50 - 60 倍!但是有点难读。
我发现，在给定 csr_matrix 格式的情况下，您可以在这里使用 indices 和 indptr 属性来发挥您的优势。
def final_app():
idx = sp_arr.indptr
arr_c = arr[:, sp_arr.indices]
out = [arr_c[:, i:j] for i, j in zip(idx[:-1], idx[1:])]
return out
最终性能在统计上是相同的(在稀疏矩阵 276538 x 33114 上改进不到 50ms)，但感觉更容易理解。 更重要的是这种方法确实包括了根本没有值的行，而以前的方法则没有。这可能看起来不重要，但对于我的用例来说它非常关键。
更新 2
回应@EelcoHoogendoorn。问题是我试图在 python 中实现的带有正则化方法的替代最小二乘并行实现的一部分。这来自经常被引用的论文 Large-scale Parallel Collaborative Filtering for the Netflix Prize这样做的正常方法是跨进程分发评分、用户和项目矩阵的副本。我认为如果我们预先构建所有项目子矩阵并将它们分发给流程，看看会发生什么可能会很有趣。这样，流程只需要分别返回一个用户或一个项目的特征列，这些列可用于更新用户和项目矩阵。
上面的问题其实是我目前实现的瓶颈。根据您的评论，在这种情况下，我认为转置并不重要，因为算法的一部分采用每个子矩阵的点积及其转置。
最佳答案
嗯，有两个选项 - np.split 或 loop comprehension。根据我的经验，我发现后者更快。但是，优先事项必须是通过尽可能多地进行预处理来在循环理解内做最少的工作。
方法 #1: 第一种方法使用 np.split -
# Get row, col indices
r,c = sp_arr.nonzero()
# Get intervaled indices for row indices.
# We need to use these to cut the column indexed input array.
idx = np.flatnonzero(r[1:] > r[:-1])+1
out = np.split(arr[:,c], idx, axis=1)
示例输出-
In [56]: [i.tolist() for i in out]
Out[56]:
[[[0.47027789, 0.3618022], [0.80198907, 0.094961]],
[[0.01321617], [0.10254934]],
[[0.82510323, 0.3618022], [0.36350616, 0.094961]],
[[0.66264852, 0.3618022], [0.65209401, 0.094961]],
[[0.66264852], [0.65209401]]]
方法 #2: 第二种应该在性能上更好，我们将重新使用以前方法中的 r,c,idx -
idx0 = np.concatenate(( [0] , idx, [r.size] ))
arr_c = arr[:,c]
out = [arr_c[:,i:j] for i,j in zip(idx0[:-1], idx0[1:])]
看，loop-comprehension 只是对已经索引的数组 arr_c 的数组进行切片。这是一个人所能达到的最低限度，因此应该很好。
运行时测试
方法-
def org_app(sp_arr, arr):
output = np.empty(sp_arr.shape[0], dtype=object)
for row in range(sp_arr.shape[0]):
output[row] = arr[:, sp_arr[row].indices]
return output
def app1(sp_arr, arr):
r,c = sp_arr.nonzero()
idx = np.flatnonzero(r[1:] > r[:-1])+1
return np.split(arr[:,c], idx, axis=1)
def app2(sp_arr, arr):
r,c = sp_arr.nonzero()
idx = np.flatnonzero(r[1:] > r[:-1])+1
idx0 = np.concatenate(( [0] , idx, [r.size] ))
arr_c = arr[:,c]
return [arr_c[:,i:j] for i,j in zip(idx0[:-1], idx0[1:])]
时间 -
In [146]: sp_arr = csr_matrix((np.random.rand(100000,100)>0.8).astype(int))
...: arr = np.random.rand(10,sp_arr.shape[1])
...:
In [147]: %timeit org_app(sp_arr, arr)
...: %timeit app1(sp_arr, arr)
...: %timeit app2(sp_arr, arr)
...:
1 loops, best of 3: 5.66 s per loop
10 loops, best of 3: 146 ms per loop
10 loops, best of 3: 105 ms per loop
关于python - 从稀疏数组的每一行快速拆分列索引数组，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/44053388/