如何计算scipy稀疏矩阵行列式而不将其变为密集？

Question

我试图找出在python中找到稀疏对称和实矩阵行列式的最快方法.使用scipy sparse模块,但真的很惊讶没有行列式功能.我知道我可以使用LU分解来计算行列式,但是没有看到一个简单的方法来执行它,因为返回scipy.sparse.linalg.splu是一个对象并且实例化一个密集的L和U矩阵是不值得的 - 我不妨sp.linalg.det(A.todense())在哪里做A我的scipy稀疏矩阵.

我也有点惊讶为什么其他人没有面对scipy中有效决定因素计算的问题.如何使用splu计算行列式？

我看着pySparse和scikits.sparse.chlmod.后者对我来说现在不实用 - 需要软件包安装,也不确定在我遇到麻烦之前代码的速度有多快.有解决方案吗 提前致谢.

Answer 1

以下是我在此处作为答案的一部分提供的一些参考资料。我认为它们解决了您试图解决的实际问题：

• Shogun 库中的实现说明
• Erlend Aune、Daniel P. Simpson：高维高斯分布中的参数估计，特别是第 2.1 节 ( arxiv:1105.5256 )
• Ilse CF Ipsen、Dean J. Lee：行列式近似( arxiv:1105.0437 )
• Arnold Reusken：大型稀疏对称正定矩阵行列式的近似( arxiv:hep-lat/0008007 )

引用将军笔记：

计算似然表达式中对数行列式项的常用技术依赖于矩阵的 Cholesky 分解，即 ?=LLT，（L 是下三角 Cholesky 因子），然后使用因子的对角线条目来计算 log(det (?))=2?ni=1log(Lii)。然而，对于稀疏矩阵，正如协方差矩阵通常那样，Cholesky 因子经常受到填充现象的影响——结果证明它们本身并不那么稀疏。因此，对于大尺寸，这种技术变得不可行，因为需要大量内存来存储所有这些不相关的非对角系数系数。虽然已经开发出排序技术来预先排列行和列以减少填充，例如近似最小度 (AMD) 重新排序，

最近的研究表明，使用复数分析、数值线性代数和贪婪图着色中的许多技术，我们可以将对数行列式近似到任意精度[Aune et. 等，2012]。主要技巧在于我们可以将 log(det(?)) 写成 trace(log(?))，其中 log(?) 是矩阵对数。

Answer 2

您可以使用scipy.sparse.linalg.splu获取分解的下 ( L) 和上 ( U) 三角矩阵的稀疏矩阵M=LU：

from scipy.sparse.linalg import splu

lu = splu(M)

行列式det(M)可以表示为：

det(M) = det(LU) = det(L)det(U)

三角矩阵的行列式只是对角项的乘积：

diagL = lu.L.diagonal()
diagU = lu.U.diagonal()
d = diagL.prod()*diagU.prod()

然而，对于大矩阵，通常会发生下溢或溢出，这可以通过使用对数来避免。

diagL = diagL.astype(np.complex128)
diagU = diagU.astype(np.complex128)
logdet = np.log(diagL).sum() + np.log(diagU).sum()

请注意，我调用复杂的算术来计算对角线上可能出现的负数。现在，logdet您可以从中恢复行列式：

det = np.exp(logdet) # usually underflows/overflows for large matrices

而行列式的符号可以直接从diagL和计算diagU（例如在实现 Crisfield 的弧长方法时很重要）：

sign = swap_sign*np.sign(diagL).prod()*np.sign(diagU).prod()

其中swap_sign是考虑 LU 分解中排列数的项。感谢@Luiz Felippe Rodrigues，可以计算出：

swap_sign = (-1)**minimumSwaps(lu.perm_r)

def minimumSwaps(arr): 
    """
    Minimum number of swaps needed to order a
    permutation array
    """
    # from https://www.thepoorcoder.com/hackerrank-minimum-swaps-2-solution/
    a = dict(enumerate(arr))
    b = {v:k for k,v in a.items()}
    count = 0
    for i in a:
        x = a[i]
        if x!=i:
            y = b[i]
            a[y] = x
            b[x] = y
            count+=1
    return count

Answer 3

解决这个问题的"标准"方法是使用cholesky分解,但是如果你不能使用任何新的编译代码,那么你就不走运了.最好的稀疏cholesky实现是蒂姆戴维斯的CHOLMOD,它是根据LGPL许可的,因此不适用于scipy本身(scipy是BSD).