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

Question

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

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

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

我正在尝试找出在运行时解决系统问题的最快方法；我不介意花时间进行预处理（例如计算 L 的因式分解）。

我最初的想法是预先计算 L 的 Cholesky 分解，然后在运行时使用 D 中的值更新分解（使用 cholupdate 进行 rank-1 更新），并通过反向替换快速解决问题。不幸的是，Cholesky 分解不像原始 L 矩阵那样稀疏，仅从磁盘加载它就需要 5.48 秒；作为对比，直接用反斜杠求解系统需要8.30s。

鉴于我的矩阵的形状，是否有任何其他方法可以建议您在运行时加快求解速度，无论预处理时间需要多长时间？

Answer 1

假设您正在处理网格（因为您提到了图像 - 尽管这不能保证），您对速度更感兴趣而不是精度（因为 5 秒对于 100 万个未知数来说似乎已经太慢了），我看到了几个选项.

首先，忘记确切的方法，例如 Cholesky (+reordering)。即使他们允许存储分解并将其重用于多个 rhs，您也可能需要存储在您的情况下似乎难以处理的巨大矩阵（我希望您使用反向 Cuthill McKee 或任何东西重新排序行/列否则-这会使分解变得很稀疏）。

根据您的边界条件，我将首先尝试使用 FFT 解决泊松问题的 Matlab poisolv，如果您想要 Dirichlet 边界条件而不是周期性边界条件，则可能进行重投影。它非常快，但可能不适合您的问题（您提到拉普拉斯矩阵+恒等式有 25 nnz：为什么？它是高阶拉普拉斯矩阵，在这种情况下，您可能对精度比我假设的更感兴趣？或者它实际上是一个与你描述的不同的问题？）。

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

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

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

当然，您可以预先计算、存储和保持分解的直接求解器也很有意义（尽管如果您的矩阵是常数，我不理解您关于 rank-1 更新的论点），因为只有反向替换仍然存在在运行时完成。但鉴于这忽略了问题的结构（规则网格，可能对有限精度结果感兴趣等），我会选择为这些情况设计的方法，例如类似傅里叶的方法或多重网格。这两种方法都可以在 GPU 上实现以获得更快的结果（回想一下，GPU 是为处理图像/纹理而量身定制的！）。

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