计算大型稀疏矩阵的幂和

Question

给定一个查询向量（one-hot-vector）q，大小为50000x1，一个大稀疏矩阵A，大小为50000 x 50000，nnz 为A，为0.3十亿，我想要计算r=(A + A^2 + ... + A^S)q（通常是4 <= S <=6）。

我可以使用循环迭代上面的方程

r = np.zeros((50000,1))
for i in range(S):
   q = A.dot(q)
   r += q

但我想要更快的方法。

首先想到的是A 可以是对称的，所以eigen decomposition would help for compute power of A。但由于A 是大型稀疏矩阵，分解后的密集矩阵与A 大小相同，这会导致性能下降（在内存和速度方面）。

还考虑了低秩近似。但是A 又大又稀疏，所以不确定r 哪个等级合适。

预先计算一些东西是完全可以的，比如A + A^2 + ... + A^S = B。但我希望最后的计算速度很快：计算 Bq 少于 40 毫秒。

是否有任何参考或论文或技巧？

Answer 1

即使矩阵不是稀疏的，迭代方法也是可行的方法。

A.dot(q) 的乘法复杂度为 O(N^2)，而 A.dot(A^i) 的计算复杂度为 O(N^3)。

q 是稀疏的（实际上比 A 稀疏得多）这一事实可能会有所帮助。 对于第一次迭代，A*q 可以计算为A[q_hot_index,:].T。

对于第二次迭代，A @ q 的预期密度与 A 的预期密度相同（大约 10%），因此稀疏地执行它仍然是好的。

对于之后的第三次迭代，A^i @ q 将是密集的。

由于您正在累积结果，因此您的 r 不是稀疏的很好，它可以防止索引操作。 有several 不同的方式来存储稀疏矩阵。我自己不能说我深入了解所有这些，但我认为csr_matrix、csc_matrix 对于通用稀疏矩阵来说是最紧凑的。

当您需要计算 P(A) 时，本征分解是好的，要计算 P(A)*q，本征分解只有在 P(A) 具有 A 大小的阶数时才会变得有利。特征分解的复杂度为O(N^3)，矩阵向量乘积的复杂度为O(N^2)，使用特征分解对度为D的多项式P(A)的求值可以在O(N^3 + N*D)中实现。

编辑：回答有关 cmets 的问题

1. “它可以防止索引操作”

假设你有一个稀疏矩阵[0,0,0,2,0,7,0]。这可以描述为((3,2), (5,7))。现在假设你将 1 分配给一个元素，它变成了[0,0,0,2,1,7,0]，它现在表示为((3,2), (4,1), (5,7))。赋值是通过插入数组来执行的，插入数组的复杂度为O(nnz)，其中nnz 是非零元素的数量。如果你有一个密集矩阵，你总是可以修改一个复杂的元素O(1)。

1. 复杂度中的 N 是多少？

它是矩阵A的行数或列数

1. 关于特征分解，你想说值得吗 计算 r 可以在 O(N^3 +N*D) 而不是 O(N^3 + N^2) 中实现

计算P(A) 将具有复杂性O(N^3 * D)（使用不同的常数），对于大矩阵，使用特征分解计算P(A) 可能是最有效的。但是P(A)x 具有O(N^2 * D) 的复杂性，因此使用特征分解计算P(A)x 可能不是一个好主意，除非你有很大的D (>N)，当涉及速度时。