稀疏 Hermitian 矩阵的 Python LDLT 分解，有吗？答案

【问题标题】：A Python LDLT factorization for sparse Hermitian matrices, is there one?稀疏 Hermitian 矩阵的 Python LDLT 分解，有吗？
【发布时间】：2019-07-08 15:27:33
【问题描述】：

我有一些大型（N×N，N=10,000 到 N=36,000,000）稀疏、复杂、厄密矩阵，通常是非奇异的，我有一个频谱切片问题。具体来说，我需要知道正特征值的确切数量。

我需要一个稀疏 LDLT 分解——有吗？理想情况下，它将是一个多前向算法，并且可以很好地并行化，并且可以选择仅计算 D 而不是上三角矩阵或置换矩阵。

我目前在 Matlab 中使用 ldl() 。这仅适用于实矩阵，所以我需要创建一个更大的实矩阵。此外，它总是计算 L 和 D。我需要一个更好的算法来适应 64GB RAM。我希望 Python 将更加可定制。（如果是这样，我将学习 Python。）我应该补充一点：每个节点可以获得 64GB RAM，并且可以获得 6 个节点。即使是一台有 64GB RAM 的机器，我也不想浪费 RAM 存储 L 只是为了删除它。

也许有人为 MUMPS（MUltifrontal Massively Parallel Solver）编写了 Python 前端？

我会使用非并行 Python 版本的 LDLT，因为我的很多研究都涉及许多中等大小的矩阵。

【问题讨论】：

+1 用于非常具体的精心制定的要求组合。所有的手指都交叉。你的问题是一场战斗，发生在计算复杂性的地狱角落。
我确实在 Python 中找到了用于稀疏矩阵的 LDLT 版本，但它不是真正的 LDLT，因为它需要半正定。也就是说，对我来说，所有特征值都是正的（或者在退化的情况下为零）的微不足道的情况。我的矩阵在正特征值和负特征值之间几乎均匀分布。

标签： python matrix parallel-processing sparse-matrix

【解决方案1】：

我需要一个更好的算法来适应 64GB 的 RAM。我希望 Python 能够更加可定制。（如果是这样，我会学习 Python。）

如果这是可能的话：

|>>> ( 2. * 8 ) * 10E3 ** 2 / 1E12            # a 64GB RAM can store matrices ~1.6 GB
0.0016                                        #        the [10k,10k] of complex64
|                                             #        or  [20k,20k] of    real64
|>>> ( 2. * 8 ) * 63E3 ** 2 / 1E12            # a 64GB RAM can store matrices in-RAM
0.063504                                      #        max [63k,63k] of type complex
|                                             #        but
|>>> ( 2. * 8 ) * 36E6 ** 2 / 1E12            #        the [36M,36M] one has
20736.0                                       # ~ 21PB of data
+0:00:09.876642                               #        Houston, we have a problem
#---------------------------------------------#--------and [2M7,2M7] has taken    
                                                                   ~ month
                                                                  on HPC cluster

研究需求很明确，但没有这样的语言（无论是 Matlab、Python、汇编程序、Julia 还是 LISP）可以将 21 PB 的数据存储到仅 64 GB 的空间中物理 RAM 存储，以尽可能快地进行复杂矩阵（具有这样的给定比例）特征值计算。我的意思是，将数据从 RAM 计算“卸载”到任何形式的 RAM 外存储是非常昂贵的（大约慢了 ~ +1E2 ~ +1E5 ~ 个数量级），任何这样的计算过程将产生仅用于第一次“读取”21 PB 元素数据的年龄。

如果您的研究有资金或赞助使用相当特定的计算设备基础设施，可能会有一些技巧来处理这些数字堆，但不要期望将 21 PB 的蠕虫（数据）放入 64 GB 的大容量（嗯，相当小）可以“免费”。

您可能会因为许多其他原因和/或动机而喜欢 Python，但不是因为任何更便宜但更快的 HPC 级并行计算，也不是因为在 64GB 设备内轻松处理 21PB 数据，也不是因为任何类型的毁灭稀疏矩阵操作的主要和巨大的[TIME]-domain 附加成本可见但在其用于计算期间。使一些 xTB 稀疏矩阵处理可行以在不到 1E2 [min] 而不是 2E3 内产生结果，我敢说我知道同时增加 [PSPACE] 数据缩放是多么困难并同时缩短经常 [EXPTIME] 处理持续时间...真正的计算复杂性的地狱角 ...稀疏矩阵表示与帮助至少享受一些潜在的[PSPACE]-储蓄相比，通常会带来更多的麻烦（同样，[SPACE] 和更多，好吧，[TIME] 更糟糕，因为新类型的处罚）

鉴于参数的范围，我可以肯定地打赌，即使是算法部分也无济于事，甚至量子计算设备的承诺将在我们的生命周期内保持不变，无法将如此巨大的参数空间扩展到 QC 退火炉基于非结构化量子驱动最小化器（处理），用于将任何合理长（短持续时间）的参数块序列转换为（有限物理尺寸）量子位场问题增强过程，正如 QC 社区目前使用的那样，谢谢LLNL 等人的研究创新。

抱歉，附近似乎没有这种魔法。

为 MUMPS 使用可用的 python 前端不会改变游戏的 HPC 问题，但如果你想使用它，是的，有几个可用的。

大规模高效的 HPC 级数字运算仍然是 [处理时间] x [（无论如何）数据表示的有效存储和检索] 的乘积问题的根本原因。

希望您将获得并享受舒适（pythonic 用户热衷于留在其中）和您希望拥有的 HPC 级性能（任何类型的后端）的正确组合。

【讨论】：

我编辑了这个问题以强调我正在使用稀疏矩阵。抱歉，如果不清楚。我可以使用 N=2,7000,000 进行 LDLT 分解，对于较大的矩阵，只能使用矩阵向量乘法估计一些内部特征值。计算在一个集群上花了一个月的时间。我正在寻找一种在不存储 L 的情况下对稀疏矩阵进行 LDLT 的方法。此外，还有一个用 Python 编写的凝聚态物理包，所以我使用 Phython 意味着我可以更轻松地与物理合作者合作。
@TerryLoring 正好相反，从一开始就提到了稀疏矩阵。您的问题表述清楚而合理，这是我的错，如果发布的答案没有给人留下正确的印象，即期望一些 [PSPACE] 节省（来自任何较少矩阵数据元素的最可行的稀疏表示），但忘记了如果试图处理来自这样一个[SPACE]-motivated 的数据，将不得不支付巨大的 [EXPTIME]来自任何快速数据存储的表示。过失，过失
这是一个微妙的问题，因为我不想计算任何特征值。我想知道有多少是阳性的。西尔维斯特惯性定律说，我们进行 LDLT 分解并计算对角线 D 的特征值。这种线性代数问题在时间上很好地适应了稀疏矩阵。我同意，以 10^-{12} 的精度计算内部特征值的扩展性很差。
如本答案所述，MUMPS 有 python 前端。我找到了 PyMUMPS，所以会尝试一下。