从稀疏数组的每一行快速拆分列索引数组
Fast splitting of array on column indices from each row of sparse array
假设我有一个稀疏数组和一个密集数组，它们的列数相同但行数较少：
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
Netflix Prize 这样做的正常方法是跨进程分发评级、用户和项目矩阵的副本。我认为如果我们预先构建所有项目子矩阵并将它们分发给流程，看看会发生什么可能会很有趣。这样，流程只需要 return 一个用户或一个项目的特征列，这些特征列可用于更新用户和项目矩阵。
上面的问题其实是我目前实现的瓶颈。根据您的评论，在这种情况下，我认为转置并不重要，因为算法的一部分采用每个子矩阵的点积及其转置。
嗯，有两个选项 - 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
另一种方法是使用groupby。
from itertools import groupby
rows, cols = sp_arr.nonzero()
out = [arr[:, [g[1] for g in group]] for _, group in groupby(zip(rows, cols), lambda x: x[0])]
numpy_indexed 包（免责声明：我是它的作者），允许你使用简单的一行代码来做这样的事情：
import numpy_indexed as npi
r, c = sp_arr.nonzero()
s = group_by(r).split(arr.T[c])
性能应该与当前接受的答案相似，但速度稍慢。请注意，s 的元素现在相对于您的原始布局转置，因此如果这是硬接口要求，则需要另一遍。
除了更加简洁和可测试之外，我提出这种方法的原因是，如果您在稍微更高的层次而不是子问题上描述您的问题，很可能仍然可以找到更优雅的解决方案你给了我们。以我的经验，拆分操作很少是任何事情的最终结果，因此如果需要后续操作，这可能会揭示更多在更高抽象级别表达逻辑的机会。
假设我有一个稀疏数组和一个密集数组，它们的列数相同但行数较少：
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 Netflix Prize 这样做的正常方法是跨进程分发评级、用户和项目矩阵的副本。我认为如果我们预先构建所有项目子矩阵并将它们分发给流程，看看会发生什么可能会很有趣。这样，流程只需要 return 一个用户或一个项目的特征列，这些特征列可用于更新用户和项目矩阵。
上面的问题其实是我目前实现的瓶颈。根据您的评论，在这种情况下，我认为转置并不重要，因为算法的一部分采用每个子矩阵的点积及其转置。
嗯，有两个选项 - 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
另一种方法是使用groupby。
from itertools import groupby
rows, cols = sp_arr.nonzero()
out = [arr[:, [g[1] for g in group]] for _, group in groupby(zip(rows, cols), lambda x: x[0])]
numpy_indexed 包（免责声明：我是它的作者），允许你使用简单的一行代码来做这样的事情：
import numpy_indexed as npi
r, c = sp_arr.nonzero()
s = group_by(r).split(arr.T[c])
性能应该与当前接受的答案相似，但速度稍慢。请注意，s 的元素现在相对于您的原始布局转置，因此如果这是硬接口要求，则需要另一遍。
除了更加简洁和可测试之外，我提出这种方法的原因是，如果您在稍微更高的层次而不是子问题上描述您的问题，很可能仍然可以找到更优雅的解决方案你给了我们。以我的经验，拆分操作很少是任何事情的最终结果，因此如果需要后续操作，这可能会揭示更多在更高抽象级别表达逻辑的机会。