如何有效地将稀疏矩阵列添加到另一个稀疏矩阵的每一列？
How would one go about adding a sparse matrix column to every column in another sparse matrix efficiently?
感谢您花时间看我的问题。
我正在尝试创建一个程序，可以将值数组添加到稀疏矩阵中的每个“列”，如下所示：
A = [[1,1,0],
[0,0,1],
[1,0,1]]
B = [[1],
[0],
[1]]
A + B = [[2,2,1],
[0,0,1],
[2,1,2]]
用稀疏矩阵的坐标格式表示如下：
A = [[0,0,1],
[2,0,1],
[0,1,1],
[1,2,1],
[2,2,1]]
B = [[0,0,1],
[2,0,1]]
A + B = [[0,0,2],
[2,0,2],
[0,1,2],
[2,1,1],
[0,2,1],
[1,2,1],
[2,2,2]]
我正在处理大型矩阵，由于其大小，它们必须以稀疏方式表示。我需要能够向矩阵中的每一列添加一列值，但使用专门处理稀疏三元组的算法。
我花了一整天的时间，实际上是连续 10 个小时，我真的很震惊，我在任何地方都找不到一个好的答案。用乘法执行这个操作简单且高效，但似乎没有任何时间和 space 内置到 scipy 或 numpy 中的高效解决方案来执行此操作，或者如果有，它将当我发现时杀了我。我试图实施一个解决方案，但最终效率极低。
基本上，我的解决方案在技术上可行，但时间效率很低，遵循以下步骤：
- 检查 A 和 B 之间行上的共享值，将相关的三元组值加在一起
- 添加 A 中的唯一值
- 在矩阵的列中为 i 添加 B_row_i、x_i、B_value_i，检查我们是否没有从 A 三元组中添加逐字值。
至少，我认为这就是它的作用...我现在完全精疲力尽，开始在编码时逐步退出。如果有人可以提出建议和快速解决方案，我们将不胜感激！
from scipy.sparse import coo_matrix
from tqdm import tqdm
class SparseCoordinates:
def __init__(self,coo_a,coo_b):
self.shape = coo_a.shape
self.coo_a_rows = coo_a.row
self.coo_a_cols = coo_a.col
self.coo_a_data = coo_a.data
self.coo_b_rows = coo_b.row
self.coo_b_cols = coo_b.col
self.coo_b_data = coo_b.data
self.coo_c_rows = []
self.coo_c_cols = []
self.coo_c_data = []
def __check(self,a,b,c,lr,lc,lv):
for i in range(len(lr)):
if lr[i] == a and lc[i] == b and lv[i] == c:
return True
return False
def __check_shared_rows(self):
for i in tqdm(range(len(self.coo_a_rows))):
for j in range(len(self.coo_b_rows)):
if self.coo_a_rows[i] == self.coo_b_rows[j]:
self.coo_c_rows.append(self.coo_a_rows[i])
self.coo_c_cols.append(self.coo_a_cols[i])
self.coo_c_data.append(self.coo_a_data[i] + self.coo_b_data[j])
def __add_unique_from_a(self):
a_unique = set(self.coo_a_rows) - set(self.coo_b_rows)
for i in tqdm(range(len(self.coo_a_rows))):
if self.coo_a_rows[i] in a_unique:
self.coo_c_rows.append(self.coo_a_rows[i])
self.coo_c_cols.append(self.coo_a_cols[i])
self.coo_c_data.append(self.coo_a_data[i])
def __add_all_remaining_from_b(self):
for i in tqdm(range(len(self.coo_b_rows))):
for j in range(self.shape[1]):
if not self.__check(self.coo_b_rows[i],j,self.coo_b_data[i],
self.coo_a_rows,self.coo_a_cols,self.coo_a_data):
self.coo_c_rows.append(self.coo_b_rows[i])
self.coo_c_cols.append(j)
self.coo_c_data.append(self.coo_b_data[i])
def add(self):
self.__check_shared_rows()
self.__add_unique_from_a()
self.__add_all_remaining_from_b()
return coo_matrix((self.coo_c_data,(self.coo_c_rows,self.coo_c_cols)),shape=self.shape)
对于 numpy 数组，由于 broadcasting，A+B 可以完成工作。 broadcasting 在核心迭代级别实现，利用了 strides。 scipy.sparse 未实现 broadcasting。如果 B 扩展为 (3,3) 矩阵以匹配 A 加法有效
In [76]: A = np.array([[1,1,0],
...: [0,0,1],
...: [1,0,1]])
...:
...: B = np.array([[1],
...: [0],
...: [1]])
In [77]: B.shape
Out[77]: (3, 1)
In [78]: A+B
Out[78]:
array([[2, 2, 1],
[0, 0, 1],
[2, 1, 2]])
稀疏：
In [79]: from scipy import sparse
In [81]: M=sparse.csr_matrix(A);N=sparse.csr_matrix(B)
矩阵乘法针对稀疏矩阵得到了很好的发展：
In [82]: M@N
Out[82]:
<3x1 sparse matrix of type '<class 'numpy.int64'>'
with 3 stored elements in Compressed Sparse Row format>
In [84]: N.T@M
Out[84]:
<1x3 sparse matrix of type '<class 'numpy.int64'>'
with 3 stored elements in Compressed Sparse Column format>
矩阵乘法用于行或列索引以及求和。
定义助手：
In [86]: O=sparse.csr_matrix([1,1,1])
In [87]: O
Out[87]:
<1x3 sparse matrix of type '<class 'numpy.int64'>'
with 3 stored elements in Compressed Sparse Row format>
In [88]: N@O
Out[88]:
<3x3 sparse matrix of type '<class 'numpy.int64'>'
with 6 stored elements in Compressed Sparse Row format>
在总和中使用：
In [89]: M+N@O
Out[89]:
<3x3 sparse matrix of type '<class 'numpy.int64'>'
with 7 stored elements in Compressed Sparse Row format>
In [90]: _.A
Out[90]:
array([[2, 2, 1],
[0, 0, 1],
[2, 1, 2]])
这个矩阵不如 M 稀疏 - 更少的 0。这会降低使用稀疏矩阵的好处。
coo_matrix 格式用于创建矩阵和构建新矩阵，与 sparse.bmat、sparse.hstack 和 sparse.vstack 一样，但 csr_matrix 是使用最多的数学。有时可以通过 intptr 数组迭代 'rows' 来进行自定义数学运算。有各种 SO 答案可以做到这一点。使用 coo 属性的数学通常不是一个好主意。
但是，因为 coo 重复项在转换为 csr 时会被求和，我可能可以对 M 和 N 的 coo 格式设置一个合理的求和。
In [91]: Mo=M.tocoo(); No=N.tocoo()
In [95]: Oo=O.tocoo()
In [98]: rows = np.concatenate((Mo.row, Oo.row))
In [99]: cols = np.concatenate((Mo.col, Oo.col))
In [100]: data = np.concatenate((Mo.data, Oo.data))
In [101]: sparse.coo_matrix((data,(rows,cols)),M.shape)
Out[101]:
<3x3 sparse matrix of type '<class 'numpy.int64'>'
with 8 stored elements in COOrdinate format>
In [102]: _.A
Out[102]:
array([[2, 2, 1],
[0, 0, 1],
[1, 0, 1]])
这里我加入了Mo和Oo的属性。我可以对 No 的属性执行相同的操作，但由于我已经拥有 Oo，我认为这是值得的。我会把它留给你。
In [103]: Oo.row,Oo.col,Oo.data
Out[103]:
(array([0, 0, 0], dtype=int32),
array([0, 1, 2], dtype=int32),
array([1, 1, 1]))
In [104]: No.row,No.col,No.data
Out[104]: (array([0, 2], dtype=int32), array([0, 0], dtype=int32), array([1, 1]))
我不知道这种 coo 方法是否值得付出努力。
由于您已经在使用 scipy 稀疏数组和 numpy，因此您可以非常有效地执行此操作。在您熟悉稀疏数组的索引之前，高效稀疏操作的实际实现可能看起来有点复杂，但是 numpy 有一些很好的工具。
您只需几行即可完成此操作，无论是 COO、CSR 还是 CSC 格式。作为快速概览，压缩的稀疏矩阵格式有 3 个一维 numpy 数组，它们实际存储数据。第一个 m.indptr 为每一行（对于 CSR）或列（对于 CSC）加上一个条目。 m.indptr[i] 给出了第 i 行开始的其他两个数组的索引，因此 s = m.indptr[i]:m.indptr[i+1] 给出了该行中所有这些索引的一部分。下一个数组 m.indices 给出未压缩轴中每个非零值的索引——CSR 的列、CSC 的行。 m.data 给出了这些位置的实际值。您可以经常利用这种格式来大大简化代码并节省时间。例如，如果您想将第 i 列中的每个非零值加 1：
- 转换为 CSC 格式
- 使用 indptr 数组查找第 i 列中非零值的索引
- 用这些索引切入数据数组，并添加到位（即 += 1）
m.data[m.indptr[i]:m.indptr[i+1]] += 1
就是这样！现在，我们可能需要为每一行或每一列添加不同的值。在 for 循环中执行该步骤对于大问题是不切实际的。一个有用的通用技术是构造一个与最初用零填充的数据数组大小相同的数组，然后使用切片和 numpy 函数（如 take、diff 和 repeat）的某种组合来准备使用 indptr 和 indices 数组的更改，然后最后一次更新数据数组。
当您向矩阵的每一列添加一个向量时，对值的修改仅取决于行索引 - 因此，如果我们能够想出一种方法来一次有效地将正确的值添加到每一行，我们很好。如果我们有一个 CSR 矩阵，我们就会确切地知道每一行中有多少个非零元素，并且它们将在数据数组中一起排序。使用 np.diff(m.indptr) 非常方便地给出每行中非零元素的数量，然后 m.data += np.repeat(vec, np.diff(m.indptr)) 给出一个与数据数组大小相同的数组， vec[0] 对第一行中的每个元素重复一次, vec[1] 对第二行中的每个元素重复一次，依此类推，然后将结果添加到数据数组中。这只需要O(nnz)的时间和内存，这是最优的，并且可以充分利用numpy的向量化和并行性而不需要额外的努力。
如果由于某些其他约束而需要根据未压缩轴中的索引更改值，则效率稍低的解决方案遵循完全相同的想法，但不是使用 repeat 和 indptr 数组，而是使用 np.take(vec, m.indices).这将创建一个数组，其元素 i 取自 m.indices[i] 给出的 vec 的索引 - 换句话说，它采用对应于每个非零元素的 row/column 索引的 vec 元素。这涉及以一种不太连续的方式访问内存以从 vec 中提取正确的值，但在实践中，由于与数据数组相比 vec 通常较小，并且数据更新是连续的，因此速度非常快。
编辑添加：sparsefuncs module in scikit-learn has some very instructive code, along with the similar sparsefuncs_fast Cython 模块。您可能会注意到我介绍的方法与行和列缩放函数中使用的代码非常相似——这就是我从中学到的。
感谢您花时间看我的问题。
我正在尝试创建一个程序，可以将值数组添加到稀疏矩阵中的每个“列”，如下所示：
A = [[1,1,0],
[0,0,1],
[1,0,1]]
B = [[1],
[0],
[1]]
A + B = [[2,2,1],
[0,0,1],
[2,1,2]]
用稀疏矩阵的坐标格式表示如下：
A = [[0,0,1],
[2,0,1],
[0,1,1],
[1,2,1],
[2,2,1]]
B = [[0,0,1],
[2,0,1]]
A + B = [[0,0,2],
[2,0,2],
[0,1,2],
[2,1,1],
[0,2,1],
[1,2,1],
[2,2,2]]
我正在处理大型矩阵，由于其大小，它们必须以稀疏方式表示。我需要能够向矩阵中的每一列添加一列值，但使用专门处理稀疏三元组的算法。
我花了一整天的时间，实际上是连续 10 个小时，我真的很震惊，我在任何地方都找不到一个好的答案。用乘法执行这个操作简单且高效，但似乎没有任何时间和 space 内置到 scipy 或 numpy 中的高效解决方案来执行此操作，或者如果有，它将当我发现时杀了我。我试图实施一个解决方案，但最终效率极低。
基本上，我的解决方案在技术上可行，但时间效率很低，遵循以下步骤：
- 检查 A 和 B 之间行上的共享值，将相关的三元组值加在一起
- 添加 A 中的唯一值
- 在矩阵的列中为 i 添加 B_row_i、x_i、B_value_i，检查我们是否没有从 A 三元组中添加逐字值。
至少，我认为这就是它的作用...我现在完全精疲力尽，开始在编码时逐步退出。如果有人可以提出建议和快速解决方案，我们将不胜感激！
from scipy.sparse import coo_matrix
from tqdm import tqdm
class SparseCoordinates:
def __init__(self,coo_a,coo_b):
self.shape = coo_a.shape
self.coo_a_rows = coo_a.row
self.coo_a_cols = coo_a.col
self.coo_a_data = coo_a.data
self.coo_b_rows = coo_b.row
self.coo_b_cols = coo_b.col
self.coo_b_data = coo_b.data
self.coo_c_rows = []
self.coo_c_cols = []
self.coo_c_data = []
def __check(self,a,b,c,lr,lc,lv):
for i in range(len(lr)):
if lr[i] == a and lc[i] == b and lv[i] == c:
return True
return False
def __check_shared_rows(self):
for i in tqdm(range(len(self.coo_a_rows))):
for j in range(len(self.coo_b_rows)):
if self.coo_a_rows[i] == self.coo_b_rows[j]:
self.coo_c_rows.append(self.coo_a_rows[i])
self.coo_c_cols.append(self.coo_a_cols[i])
self.coo_c_data.append(self.coo_a_data[i] + self.coo_b_data[j])
def __add_unique_from_a(self):
a_unique = set(self.coo_a_rows) - set(self.coo_b_rows)
for i in tqdm(range(len(self.coo_a_rows))):
if self.coo_a_rows[i] in a_unique:
self.coo_c_rows.append(self.coo_a_rows[i])
self.coo_c_cols.append(self.coo_a_cols[i])
self.coo_c_data.append(self.coo_a_data[i])
def __add_all_remaining_from_b(self):
for i in tqdm(range(len(self.coo_b_rows))):
for j in range(self.shape[1]):
if not self.__check(self.coo_b_rows[i],j,self.coo_b_data[i],
self.coo_a_rows,self.coo_a_cols,self.coo_a_data):
self.coo_c_rows.append(self.coo_b_rows[i])
self.coo_c_cols.append(j)
self.coo_c_data.append(self.coo_b_data[i])
def add(self):
self.__check_shared_rows()
self.__add_unique_from_a()
self.__add_all_remaining_from_b()
return coo_matrix((self.coo_c_data,(self.coo_c_rows,self.coo_c_cols)),shape=self.shape)
对于 numpy 数组，由于 broadcasting，A+B 可以完成工作。 broadcasting 在核心迭代级别实现，利用了 strides。 scipy.sparse 未实现 broadcasting。如果 B 扩展为 (3,3) 矩阵以匹配 A 加法有效
In [76]: A = np.array([[1,1,0],
...: [0,0,1],
...: [1,0,1]])
...:
...: B = np.array([[1],
...: [0],
...: [1]])
In [77]: B.shape
Out[77]: (3, 1)
In [78]: A+B
Out[78]:
array([[2, 2, 1],
[0, 0, 1],
[2, 1, 2]])
稀疏：
In [79]: from scipy import sparse
In [81]: M=sparse.csr_matrix(A);N=sparse.csr_matrix(B)
矩阵乘法针对稀疏矩阵得到了很好的发展：
In [82]: M@N
Out[82]:
<3x1 sparse matrix of type '<class 'numpy.int64'>'
with 3 stored elements in Compressed Sparse Row format>
In [84]: N.T@M
Out[84]:
<1x3 sparse matrix of type '<class 'numpy.int64'>'
with 3 stored elements in Compressed Sparse Column format>
矩阵乘法用于行或列索引以及求和。
定义助手：
In [86]: O=sparse.csr_matrix([1,1,1])
In [87]: O
Out[87]:
<1x3 sparse matrix of type '<class 'numpy.int64'>'
with 3 stored elements in Compressed Sparse Row format>
In [88]: N@O
Out[88]:
<3x3 sparse matrix of type '<class 'numpy.int64'>'
with 6 stored elements in Compressed Sparse Row format>
在总和中使用：
In [89]: M+N@O
Out[89]:
<3x3 sparse matrix of type '<class 'numpy.int64'>'
with 7 stored elements in Compressed Sparse Row format>
In [90]: _.A
Out[90]:
array([[2, 2, 1],
[0, 0, 1],
[2, 1, 2]])
这个矩阵不如 M 稀疏 - 更少的 0。这会降低使用稀疏矩阵的好处。
coo_matrix 格式用于创建矩阵和构建新矩阵，与 sparse.bmat、sparse.hstack 和 sparse.vstack 一样，但 csr_matrix 是使用最多的数学。有时可以通过 intptr 数组迭代 'rows' 来进行自定义数学运算。有各种 SO 答案可以做到这一点。使用 coo 属性的数学通常不是一个好主意。
但是，因为 coo 重复项在转换为 csr 时会被求和，我可能可以对 M 和 N 的 coo 格式设置一个合理的求和。
In [91]: Mo=M.tocoo(); No=N.tocoo()
In [95]: Oo=O.tocoo()
In [98]: rows = np.concatenate((Mo.row, Oo.row))
In [99]: cols = np.concatenate((Mo.col, Oo.col))
In [100]: data = np.concatenate((Mo.data, Oo.data))
In [101]: sparse.coo_matrix((data,(rows,cols)),M.shape)
Out[101]:
<3x3 sparse matrix of type '<class 'numpy.int64'>'
with 8 stored elements in COOrdinate format>
In [102]: _.A
Out[102]:
array([[2, 2, 1],
[0, 0, 1],
[1, 0, 1]])
这里我加入了Mo和Oo的属性。我可以对 No 的属性执行相同的操作，但由于我已经拥有 Oo，我认为这是值得的。我会把它留给你。
In [103]: Oo.row,Oo.col,Oo.data
Out[103]:
(array([0, 0, 0], dtype=int32),
array([0, 1, 2], dtype=int32),
array([1, 1, 1]))
In [104]: No.row,No.col,No.data
Out[104]: (array([0, 2], dtype=int32), array([0, 0], dtype=int32), array([1, 1]))
我不知道这种 coo 方法是否值得付出努力。
由于您已经在使用 scipy 稀疏数组和 numpy，因此您可以非常有效地执行此操作。在您熟悉稀疏数组的索引之前，高效稀疏操作的实际实现可能看起来有点复杂，但是 numpy 有一些很好的工具。
您只需几行即可完成此操作，无论是 COO、CSR 还是 CSC 格式。作为快速概览，压缩的稀疏矩阵格式有 3 个一维 numpy 数组，它们实际存储数据。第一个 m.indptr 为每一行（对于 CSR）或列（对于 CSC）加上一个条目。 m.indptr[i] 给出了第 i 行开始的其他两个数组的索引，因此 s = m.indptr[i]:m.indptr[i+1] 给出了该行中所有这些索引的一部分。下一个数组 m.indices 给出未压缩轴中每个非零值的索引——CSR 的列、CSC 的行。 m.data 给出了这些位置的实际值。您可以经常利用这种格式来大大简化代码并节省时间。例如，如果您想将第 i 列中的每个非零值加 1：
- 转换为 CSC 格式
- 使用 indptr 数组查找第 i 列中非零值的索引
- 用这些索引切入数据数组，并添加到位（即 += 1）
m.data[m.indptr[i]:m.indptr[i+1]] += 1
就是这样！现在，我们可能需要为每一行或每一列添加不同的值。在 for 循环中执行该步骤对于大问题是不切实际的。一个有用的通用技术是构造一个与最初用零填充的数据数组大小相同的数组，然后使用切片和 numpy 函数（如 take、diff 和 repeat）的某种组合来准备使用 indptr 和 indices 数组的更改，然后最后一次更新数据数组。
当您向矩阵的每一列添加一个向量时，对值的修改仅取决于行索引 - 因此，如果我们能够想出一种方法来一次有效地将正确的值添加到每一行，我们很好。如果我们有一个 CSR 矩阵，我们就会确切地知道每一行中有多少个非零元素，并且它们将在数据数组中一起排序。使用 np.diff(m.indptr) 非常方便地给出每行中非零元素的数量，然后 m.data += np.repeat(vec, np.diff(m.indptr)) 给出一个与数据数组大小相同的数组， vec[0] 对第一行中的每个元素重复一次, vec[1] 对第二行中的每个元素重复一次，依此类推，然后将结果添加到数据数组中。这只需要O(nnz)的时间和内存，这是最优的，并且可以充分利用numpy的向量化和并行性而不需要额外的努力。
如果由于某些其他约束而需要根据未压缩轴中的索引更改值，则效率稍低的解决方案遵循完全相同的想法，但不是使用 repeat 和 indptr 数组，而是使用 np.take(vec, m.indices).这将创建一个数组，其元素 i 取自 m.indices[i] 给出的 vec 的索引 - 换句话说，它采用对应于每个非零元素的 row/column 索引的 vec 元素。这涉及以一种不太连续的方式访问内存以从 vec 中提取正确的值，但在实践中，由于与数据数组相比 vec 通常较小，并且数据更新是连续的，因此速度非常快。
编辑添加：sparsefuncs module in scikit-learn has some very instructive code, along with the similar sparsefuncs_fast Cython 模块。您可能会注意到我介绍的方法与行和列缩放函数中使用的代码非常相似——这就是我从中学到的。