快速稀疏矩阵乘法

Question

对于类，我必须为稀疏矩阵编写我自己的线性方程求解器。我可以自由地为稀疏矩阵使用任何类型的数据结构，我必须实现几个解决方案，包括共轭梯度。

我想知道是否有一种著名的方法来存储稀疏矩阵，以便与向量相乘相对较快。

现在我的稀疏矩阵基本上实现了一个包装的std::map< std::pair<int, int>, double>，它存储数据，如果有的话。这就把一个复杂度为O(n²)的矩阵的乘法变换成了O(n²log(n))因为我需要查找每个矩阵元素。我已经研究了耶鲁稀疏矩阵格式，似乎一个元素的检索也是在O(log(n))，所以我不确定它是否会快得多。

作为参考，我有一个800x800的矩阵，填充了5000个条目。用共轭梯度法求解这样的方程组大约需要450秒。

你认为有可能用另一种数据结构做得更快吗?

谢谢!

Answer 1

最常见的选择是CSC或CSR存储。这些都是矩阵-向量乘法的有效方法。如果你必须自己做的话，编写这些乘法例程也很容易。

也就是说，耶鲁存储也产生非常有效的矩阵向量乘法。如果您正在执行矩阵元素查找，那么您就误解了如何使用该格式。我建议你学习一些标准的稀疏库来学习矩阵-向量乘法是如何实现的。

即使使用你当前的存储，你也可以以0 (n)的复杂度执行矩阵乘法。我所见过的所有稀疏矩阵向量乘法算法都可以归结为相同的步骤。例如，考虑y = Ax

1. 将结果向量y归零
2. 为矩阵a的非零元素初始化迭代器。
3. 获取矩阵A[i,j]的下一个非零元素。注意，i,j的模式不遵循常规模式。它只是反映了A的非零元素的存储顺序。
4. y[i] += A[i,j]*x[j]
5. 如果有更多的A元素，转到3

我怀疑你在写经典的double for循环密集乘法代码:

for (i=0; i<N; i++)
    for (j=0; j<N; j++)
        y[i] += A[i,j]*x[j]

，这就是导致您执行查找的原因。

但我不是建议你坚持你的std::map存储。这不会是非常有效的。我推荐CSC主要是因为它使用最广泛。