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

Question

在实现图像处理算法时，我必须求解一个大型线性系统的表单A*x=b，其中：

• 矩阵A=L+D是拉普拉斯矩阵L和对角矩阵D的和。
• 拉普拉斯矩阵L是稀疏的，每行约有25个非零矩阵。
• 系统很大，与输入图像中的像素一样多未知数(通常超过100万)。

Laplacian矩阵L在算法的连续运行中不发生变化，可以在预处理过程中构造该矩阵，并可能计算其因式分解。算法每运行一次，对角线矩阵D和右侧向量b都会发生变化.

我试图找出在运行时解决系统问题的最快方法；我不介意花时间在预处理上(例如，用于计算L的因式分解)。

我最初的想法是预先计算L的Cholesky因式分解，然后在运行时用D的值更新因式分解(用cholupdate进行秩-1更新)，然后用反向替换快速解决问题。不幸的是，Cholesky分解并不像原始L矩阵那样稀疏，仅仅从磁盘加载它就需要5.48s；相比之下，用反斜杠直接求解系统需要8.30。

考虑到我的矩阵的形状，是否还有其他方法可以在运行时加速求解，不管预处理时间需要多长时间？

Answer 1

假设您正在网格上工作(因为您提到了图像--尽管这并不一定)，您更感兴趣的是速度而不是精度(因为对于100万未知数来说，5s似乎已经太慢了)，我看到了几个选项。

首先，忘记精确的方法，比如Cholesky (+重新排序)。即使它们允许存储因式分解并将其重用于多个rhs，您也可能需要存储在您的情况下似乎难以处理的巨大矩阵(我希望您正在使用反向Cuthill McKee或其他任何东西重新排序行/列)，这将在很大程度上促进分解。

根据你的边界条件，我首先尝试一个Matlab poisolv，它用快速傅立叶变换来解决泊松问题，如果你想要Dirichlet边界条件而不是周期边界条件，可以重新投影。它非常快，但可能不适合您的问题(您提到有25 nnz的拉普拉斯matrix+identity :为什么？这是一个高阶拉普拉斯矩阵吗?在这种情况下，你可能比我假设的更感兴趣的是精确性？或者--事实上，这是一个不同于你描述的问题?)

然后，您可以尝试多重网格求解，这是非常快的图像和平滑的问题。您可以对多重网格的每个迭代和每个级别使用简单的松弛方法，也可以使用更高级的方法(例如，预处理的共轭梯度par级别)。或者，您可以在不使用多重网格的情况下进行更简单的预条件共轭梯度(甚至SSOR)，如果您只对近似解感兴趣，则可以在完全收敛之前停止迭代。

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

• 如果你想要一个近似的问题，你可以在收敛之前停止。
• 您仍然可以重用其他结果来初始化您的解决方案(例如，如果您的不同运行对应于视频的不同帧，那么使用前一个帧的解决方案作为下一个帧的初始化将是有意义的)。

当然，您可以预先计算、存储和保持因式分解的直接求解器也是有意义的(虽然我不理解如果矩阵是常数的话，您的1级更新的论点)，因为在运行时只需要进行反向替换。但是，考虑到这忽略了问题的结构(规则网格，可能对有限的精度结果感兴趣等)，我会选择为这些情况设计的方法，例如Fourier类方法或多重网格。这两种方法都可以在GPU上实现，以获得更快的结果(回想一下，GPU是为处理图像/纹理而量身定做的！)

最后，您可以从scicomp.stackexchange得到更有针对性的数值分析的有趣答案。