Eigen - 直接计算巨大稀疏矩阵的对数行列式答案

【问题标题】：Eigen - directly compute log determinant of huge sparse matrixEigen - 直接计算巨大稀疏矩阵的对数行列式
【发布时间】：2016-11-20 11:49:37
【问题描述】：

我想计算一个非常大的矩阵 (5e6 x 5e6) 的对数行列式。然而，它是高度稀疏的——每行只有 6 个非零条目（计算对角线的有 7 个）。它也是对称正定的。

在 Eigen 中，我尝试使用 Cholesky 分解：SimplicialLDLT<SparseMatrix<double>>，然后对对角线的对数值求和（SimplicialLDLT::vectorD() 可访问）但是分解运行了很长时间而没有完成。有更好的方法吗？我实际上不需要任何类型的分解，只需要对数行列式本身（或一个好的估计）。

【问题讨论】：

对于这样的矩阵，SimplicialLDLT 应该非常快，大约为 1-2 秒。确保在编译器优化开启的情况下进行编译，并且正确填充了矩阵。
你能确定你没有点击交换吗？我问是因为，在研究算法之前，也许出于某种原因，Eigen 试图密集填充一个巨大的矩阵？
另外，非零条目的大小按什么顺序排列？我正在尝试创建一个类似的数组进行测试，但我不确定是否应该用统一随机数据或正常随机数据或其他填充非零条目。
所以一般来说，对称矩阵的 LU 分解会沿着 L 或 U 的对角线有正负元素。你是否注意，当你记录负数时，你要么保持绝对值吗？还是使用具有复杂功能的日志来返回一个复数，然后再求和？
@AhmedFasih 所有非零项非对角线均为-1。在对角线上，它们是 1-5 之间的小值。也许你是对的，它试图让它变得稠密，但我不知道该怎么说。

标签： linear-algebra sparse-matrix eigen

【解决方案1】：

您期望多快？您可以尝试 3.3-beta1 版本中的所有稀疏求解器，然后找出哪个更快解决您的问题。

https://eigen.tuxfamily.org/dox-devel/group__TopicSparseSystems.html

【讨论】：

【解决方案2】：

我不妨把这个作为答案，这样我就可以展示一个数字了。

首先， Eigen 的 documentation on sparse solvers 说 SimplicialLDLT 是“推荐用于非常稀疏且不太大的问题”。你的问题很稀疏，但也很大。

第二， SimplicialLDLT 需要输入不仅是对称的（“selfadjoint”）而且是正定的，你的几乎肯定不是（除非你有理由不这么想？）。

SimplicialLDLT 可能会花费大量时间计算 Cholesky 因式分解，因为您的矩阵不是正定矩阵，所以它无法成功找到。

这让我想到了第三点。我生成了in Matlab ?，您的问题的一个小版本：1e4 x 1e4 稀疏矩阵，对称，对角线上的小整数介于 1 和 5 之间，每行有六个其他 -1 条目。 Matlab 的 LDLT 等价物是 ldl 并且奇怪的是不需要正定矩阵，因此它在几秒钟内就通过了我的 1e4 x 1e4 示例（生成稀疏矩阵比将其分解为 LDL 需要更长的时间')。

这是下三角因子的稀疏模式（通过 Matlab 中的spy(L)）：

它有 1e7 个非零元素，占用 162 MB RAM。回想一下，这是一个 1e4 x 1e4 问题。如果内存使用量与矩阵的长度成线性关系（1e4 → 5e6），那么您看到的内存使用量接近 80 GB。相反，如果它与元素的数量（1e4^2 → 5e6^2）一起缩放，则需要 38 TB RAM……这个分析都不是结论性的——很可能是通过 5e6 缩放到 5e6 大大增加了 LDL 的稀疏性。因素，但这可以解释为什么 Eigen 挂起。如 cmets 中所述，检查您的交换文件是否抖动。

第四个问题是，对于我的 1e4 x 1e4 测试示例，我在下三角 L 矩阵的对角线上有很多完全为零的条目，所以整个稀疏矩阵为零到双精度，对数或无对数。

【讨论】：