4

在我的图像处理算法的实现中，我必须求解形式为 的大型线性系统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。

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

4

1 回答 1

2

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

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

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

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

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

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

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

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

于 2013-02-01T23:57:00.060 回答