c++ - C++ 中半定矩阵的 Cholesky 分解

Question

我要实现一个函数来在 C++ 中对半定矩阵进行 Cholesky 分解，我想知道是否有一个库/任何已经优化过的东西。它的工作方式与本文所述的一样或类似:

http://www.sciencedirect.com/science/article/pii/S0096300310012713

这是正定的示例，但不适用于半正定:http://en.wikipedia.org/wiki/Cholesky_decomposition#The_Cholesky.E2.80.93Banachiewicz_and_Cholesky.E2.80.93Crout_algorithms

程序必须是 C++ 语言，没有 C/FORTRAN 库，(想一想尖锐的野蛮老板给出的说明)这意味着 ATLAS、LAPACK 等。出局。我查看了 MTL + Boost，但它们只适用于正定矩阵。是否有任何我没有找到的库，甚至是已经编写的单个函数？

Answer 1

Cholesky 半定矩阵分解的问题在于 1) 它不是唯一的 2) Crout 的算法失败。

分解的存在性通常通过限制性论证非建设性地证明(如果 M_n -> M，并且 M_n = U^t_n U_n，则 ||U_n|| = ||M_n||^1/2，其中||.|| = Hilbert-Schmidt范数，U_n为有界序列。提取子序列求极限U满足U^t U = M和U三角。)

我发现，在我感兴趣的情况下，将对角线元素乘以 1 + epsilon 是令人满意的，其中 epsilon 较小(取机器 epsilon 的几千倍)以给出完全可以接受的分解。

确实，如果 M 是半正定的，那么对于每个 epsilon > 0，M + epsilon I 是正定正的。

当 epsilon 变为零时方案收敛，您可以考虑计算多个 epsilon 的分解并执行 Richardson 外推。

至于正定情况，您可以自己实现 Crout 的算法(Numerical Recipes 中有示例代码)，但我强烈建议您不要自己编写，并建议改用 LAPACK。

如果您的老板担心 LAPACK 的潜在实现不佳，这可能涉及让您的老板为英特尔 MKL 付费。大多数时候我听到这样的发言，理由是“但是我们无法控制代码，我们确实想自己写，以便在出现问题时进行调试”。愚蠢的论点。 LAPACK 已有 40 年历史并经过全面测试。

要求不使用 LAPACK 就像要求不使用正弦、余弦和对数的标准库一样愚蠢。