如何高效地求解带有拉普拉斯和对角矩阵的线性系统？

Question

如何高效地求解带有拉普拉斯和对角矩阵的线性系统？

matlabsystemsparse-matrixfactorization

4

在我的图像处理算法实现中，我需要解决一个大型线性系统的形式为 A*x=b，其中：

- 矩阵 A=L+D 由拉普拉斯矩阵 L 和对角矩阵 D 相加得到 - 拉普拉斯矩阵 L 是稀疏的，每行大约有25个非零元素 - 这个系统很大，未知量与输入图像中的像素数量相同（通常>100万）

拉普拉斯矩阵 L 在算法的连续运行之间不会改变；我可以在预处理中构建这个矩阵，并可能计算其分解。对角矩阵 D 和右侧向量 b 在每次运行算法时都会改变。

我想找出在运行时解决此问题的最快方法；我不介意在预处理上花费时间（例如计算 L 的分解）。

我的初始想法是预先计算 L 的 Cholesky 分解，然后在运行时使用来自 D 的值更新分解（使用 cholupdate 进行秩1更新），并使用回代快速解决问题。不幸的是，Cholesky 分解不像原始的 L 矩阵那样稀疏，仅从磁盘加载就需要5.48秒；相比之下，使用 backslash 直接解决该系统需要8.30秒。

鉴于我的矩阵形状，您是否推荐任何其他方法来加速运行时的求解，无论预处理时间有多长？

- Asmo

1

试试奇异值分解（SVD），虽然我不知道对于那么大量的数据它有多快... - Ander Biguri

我不明白为什么在这里应用SVD是个好主意，考虑到我的矩阵形式。 - Asmo

是的，也许你是对的，我不是专家。只是以防万一你没有考虑到那个 :) - Ander Biguri

看一下这个链接：http://scicomp.stackexchange.com/questions/1001/how-does-the-matlab-backslash-operator-solve-ax-b-for-square-matrices - johnish

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- nbonneel · Accepted Answer

假设您正在使用网格（因为您提到了图像-尽管这并不保证），而且您更关心速度而不是精度（因为5秒对于100万个未知数似乎已经太慢了），我看到有几个选项。

首先，忘记关于Cholesky（+重新排序）等精确方法。即使它们允许存储分解并将其重用于多个rhs，您也可能需要存储巨大的矩阵，在您的情况下似乎是难以处理的（我希望您使用逆Cuthill McKee或其他任何东西重新排序行/列-这样可以大大稀疏化因式分解）。

根据您的边界条件，我首先会尝试使用Matlab的poisolv，该函数使用FFT解决泊松问题，并进行可能的重投影，如果您想要的是Dirichlet边界条件而不是周期性的话。它非常快，但可能不适合您的问题（您提到Laplacian矩阵+identity具有25个非零元：为什么？这是高阶Laplace矩阵吗，在这种情况下，您可能更关心精度而不是我所假设的内容？还是实际上是与您描述的问题不同的问题？）。

然后，您可以尝试多重网格求解器，对于图像和平滑问题非常快速。您可以在每个迭代和每个多重网格级别上使用简单的松弛方法，或者使用更高级的方法（例如，每个级别的预处理共轭梯度）。另外，您可以使用较简单的预处理共轭梯度（甚至是SSOR）而不使用多重网格，如果您只对近似解感兴趣，则可以在完全收敛之前停止迭代。

我对迭代求解器的论点是：

- 如果您想要近似解，则可以在收敛之前停止。 - 您仍然可以重复使用其他结果来初始化您的解决方案（例如，如果您的不同运行对应于视频的不同帧，则使用先前帧的解决方案作为下一个帧的初始化将有一定意义）。

当然，对于可以预先计算、存储和保留分解的直接求解器也是有意义的（尽管我不理解如果矩阵是常数时进行秩-1更新的论点），因为在运行时只剩下回代需要完成。但是，考虑到这忽略了问题的结构（一个常规网格，可能对有限精度结果感兴趣等），我会选择专为这些情况设计的方法，如类傅里叶方法或多重网格方法。这两种方法都可以在GPU上实现，以获得更快的结果（请记住，GPU更适合处理图像/纹理！）。

最后，您可以从scicomp.stackexchange中获得有趣的答案，该网站更针对数值分析。