稀疏矩阵 Cholesky 分解的复杂性答案

【问题标题】：Complexity of Sparse Matrix Cholesky decomposition稀疏矩阵 Cholesky 分解的复杂性
【发布时间】：2021-09-30 05:26:47
【问题描述】：

我很难找到以下问题的直接答案：

如果计算 nxn 正定对称矩阵 A 的 Cholesky 分解，即因子 A=LL^T 且 L 为下三角矩阵，则复杂度为 O(n^3)。 对于稀疏矩阵，显然有更快的算法，但要快多少？

对于这样一个包含 m

编辑：我的矩阵也近似于主对角线（只有对角线和下方和上方的一些相邻对角线是非零的）。

P.S 我最终对 Julia 或 Python 中的实现感兴趣。 Python 有 sksparse.cholmod 模块 (https://scikit-sparse.readthedocs.io/en/latest/cholmod.html) 但我不清楚他们使用的是什么算法以及它的复杂性是什么。不确定 Julia，如果有人可以告诉我。

【问题讨论】：

维基百科有一个关于 Cholesky 分解的页面。 “当它适用时，它的效率大约是 LU 分解的两倍。” ——
代码可以在不改变复杂性的情况下更快。实际速度通常是某种固定设置代码加上取决于元素数量的每个元素复杂度的混合。但是对于稀疏，时间方式也取决于稀疏性，即相对于总数的非零元素的数量。对于像 python 这样的语言，增加了解释与编译的因素。我希望看到理论论文中讨论的复杂性具有狭义的比较；它在工作代码中的用处不大。
稀疏求解器往往不适合复杂性分析，因为它们受到稀疏结构的严重影响。不仅仅是有多少非零，而是它们在哪里。我倾向于通过经验评估性能，例如：dl.acm.org/doi/10.1145/1236463.1236465

标签： python julia sparse-matrix gaussian-process

【解决方案1】：

这只能在 P=NP 的情况下准确地回答任意矩阵......因此一般无法回答。时间复杂度取决于使用的填充减少排序，它试图获得一个 NP 难题的近似解。

但是，对于来自规则方形 2D 或 3D 网格的矩阵的非常特殊的情况，有一个答案。在这种情况下，嵌套剖析给出了渐近最优的排序。对于 2D s-by-s 网格，矩阵的维度为 n = s^2，我认为大约有 5n 个条目。在这种情况下，L 有 31*(n log2(n)/8)+O(n) 个非零，工作量是 829*(n^(3/2))/84+O(n log n)。对于 n = s^3 的 3D s-by-s-by-s 网格，L 中有 O(n^(4/3)) 个非零值，并且需要 O(n^2) 个操作来计算 L。

【讨论】：

很好的答案！您能否分享一份参考资料以进一步了解常规网格案例？
是的，见archive.siam.org/books/fa02。
还有link.springer.com/article/10.1007/BF01396660
原论文在这里：doi.org/10.1137/0710032

【解决方案2】：

Python 库 Numpy（数值 Python）还有一个用于 cholesky 的模块 - np.linalg.cholesky，我提供了文档的链接，尽管我不确定这是否能回答问题，可能需要一些实验。

【讨论】：

阿格尼，谢谢。但我最感兴趣的是复杂性。

【解决方案3】：

值得研究一个不完整的 Cholesky 分解，它有多种变体，但通常要么只计算输入中非零的三角因子中的条目，要么使用分解的低秩近似。从您的问题中不清楚您为什么对渐近复杂性感兴趣，但是使用高斯过程标签我可以猜到您正在分解协方差核矩阵并在推理过程中反复求解线性系统。在这类应用程序中，不完全分解最常用作预处理器 - 虽然它不精确，但它非常有效，易于增量更新，并且可以大大加快求解器的速度。

然而，在应用程序中，您问题的答案很可能与最适合您的目的的方法无关。对于不完全分解有 $O(n log(n)^k)$ 时间算法，由于与具有最低指数的矩阵乘法算法相同的原因，这些算法没有实际用途。知道您的应用程序究竟需要什么将告知您的选择，但您已经拥有最好的工具来找到 - 花几分钟时间编写一些代码来生成合成数据或对不同大小的真实数据进行采样并在与您对应用程序感兴趣的方式相同。如果您要进行一次分解并求解多个系统，则进行分解的时间可能会因求解而相形见绌。如果您要进行多次分解，尤其是每次都从头开始，则运行时间的常数和线性因素将产生更大的影响。另外，稀疏模式和内核本身可以对相同算法的性能产生巨大影响。构造具有完全密集的协方差矩阵、三对角精度矩阵和 Cholesky 分解的核并不困难，该分解恰好是矩阵的下三角部分，恰好可以按对角矩阵进行缩放（1d 中的指数核同时具有所有 3 个） .首先配置文件，最后优化。

【讨论】：