python - 将稀疏矩阵转换为密集矩阵并获得完整的特征值

Question

最近我正在研究一个需要解决的问题对角化一个巨大的厄尔米特矩阵以获得所有特征值。目前我正在使用 Mathematica 来完成这项工作。

但由于内存限制不适用当矩阵大小接近 (2^15,2^15) 时，对角化大约需要 32 GB 的内存。

我尝试通过从 Mathematica 导入矩阵来使用 Python，

import numpy as np
from scipy.io import mmread
from scipy.sparse import csc_matrix

#importing sparse matrix to save space
h = mmread("h.mtx")
h = csc_matrix(h)
#diagonlizing the dense one
ev = np.linalg.eigvalsh(h.todense())

它可以工作，但不幸的是比 Mathematica 慢了一个数量级。

那么，是否还有其他可能的解决方案，比如 C++？

我对 C++ 一无所知，所以我想最简单的方法可能是导入矩阵到 C++ 和对角化。

谢谢!

Answer 1

使用这个矩阵运行一些初步测试:

http://math.nist.gov/MatrixMarket/data/NEP/h2plus/qc2534.html

我确定转换为密集不会占用太多时间。特征值计算确实如此。

Numpy 使用高度优化的 Lapack 例程进行计算。这些与您在 C++ 中使用的相同。因此 C++ 不会给你太多的加速。如果您想加速使用稀疏性作为属性，请使用更好的计算机或切换到分布式矩阵存储(这里需要大量工作)。

P.S:如果您为大学项目执行此操作，您可能需要四处看看您的大学是否有某种集群。集群节点通常有很多内存。如果没有，请检查亚马逊的 AWS EC2 或谷歌的计算引擎以查找具有大量内存的实例。

编辑:

这里 Wolfram 说明了 Mathematica 在幕后所做的事情:http://reference.wolfram.com/language/tutorial/LinearAlgebraAppendix.html#83486633

Arpack 是一个 (arnoldi) 子空间求解器，只为您提供最高或最低的 k-特征值，ATLAS 只是一个 Lapack 实现，其余似乎用于求解线性系统。

所有为您提供完整特征谱的方法都需要对 NxN 矩阵进行矩阵分解。如果您只需要 k 个 vector ，可以使用一些方法将其简化为 k x k 矩阵的分解。

Arpack 有现代替代品(http://slepc.upv.es/ 或 MKL 附带的替代品)，但它们都为您提供子空间。