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

Posted 2023-03-12

【中文标题】如何计算 scipy 稀疏矩阵行列式而不将其变为密集？

【英文标题】：How to compute scipy sparse matrix determinant without turning it to dense?

【发布时间】：2013-10-07 02:34:15

【问题描述】：

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

我也有点惊讶为什么其他人没有遇到 scipy 中有效行列式计算的问题。如何使用splu 来计算行列式？

我查看了pySparse 和scikits.sparse.chlmod。后者现在对我来说不实用 - 需要安装包，并且在我陷入所有麻烦之前也不确定代码有多快。 有什么解决办法吗？提前致谢。

【问题讨论】：

您希望如何做到这一点取决于矩阵中的数据。如果您只想知道矩阵是否是奇异的，那么 eigs(A, 1, which="SM", return_eigenvectors=False) 或 svds(A, 1, which="SM", return_singular_vectors=False) 之类的东西可能是您的矩阵是否奇异的一个很好的指标。我不愿意说它总是会起作用......

首先，eigs只能返回

是的，我建议专门针对您试图查看数组是否为单数的情况。如果特征值之一为 0，则行列式也应为零。不过，这仅在数据不太疯狂的情况下才有效，因为您想查看您是否在 0 的某个容差范围内。

你说你查看了pySparce，它有一个superlu 接口，它通过部分旋转来实现LU 分解，为什么它不能满足你的需求？

因为，正如我在上面的问题中提到的，实例化 L 和 U 矩阵并不是解决问题的实用方法，L 和 U 将是密集的。稀疏矩阵方法的重点不是那样做。

【参考方案1】：

以下是我在回答 here 时提供的一些参考资料。 我认为它们解决了您要解决的实际问题：

notes 用于 Shogun 库中的实现 Erlend Aune、Daniel P. Simpson：高维高斯分布中的参数估计，特别是第 2.1 节 (arxiv:1105.5256) Ilse C.F. 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 等。等人，2012]。主要技巧在于我们可以将 log(det(Σ)) 写为 trace(log(Σ))，其中 log(Σ) 是矩阵对数。

【讨论】：

Shogun 图书馆链接已损坏。

【参考方案2】：

解决此问题的“标准”方法是使用胆汁分解法，但如果您不习惯使用任何新的编译代码，那么您就不走运了。最好的稀疏 cholesky 实现是 Tim Davis 的 CHOLMOD，它在 LGPL 下获得许可，因此在 scipy 中不可用（scipy 是 BSD）。

【讨论】：

这很简单，只需将 L 和 U 矩阵上的对角线项相乘即可。但是如果 scipy.sparse.linalg.splu 函数返回密集的 L 和 U 矩阵，那么这并不能满足 OP 的内存需求。

我知道真正的对称矩阵行列式是使用 Cholesky 计算的——事实上，我在上面的帖子中确实提到了 scikit.sparse 中的 chlmod。重点是，不要实例化密集的上下三角矩阵。我想我在重复自己。请查看我的其他 cmets。

【参考方案3】：

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

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()

但是，对于large matrices underflow or overflow commonly occurs，可以通过使用对数来避免。

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 = 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

【讨论】：

可以算作，例如：number_of_permutations = ((lu.perm_r != np.arange(lu.perm_r.size)).sum() + (lu.perm_c != np.arange(lu.perm_c.size)).sum())//2

@LuizFelippeRodrigues 好问题。也许他必须包含在建议的算法中。您是否比较了考虑和不考虑 SPLU 排列次数的结果？

它确实对我的测试产生了影响：不幸的是，您建议的最后一行似乎并没有为我带来正确的符号。但是，我之前关于如何从lu.perm_r/lu.perm_c 计算排列数的建议是错误的。我正在尝试提出解决方案。

@LuizFelippeRodrigues 这方面有什么进展吗？

如果我将代码给出的符号乘以 (-1) 到 splu 使用的行排列数，我会在测试中得到正确的结果。我已经在this gist 中实现了这一点。但是，我仍然想知道是否有更有效的解决方案..

【参考方案4】：

使用 SuperLU 和 CHOLMOD 在 N=1e6 附近的稀疏三对角线 (-1 2 -1) 行列式开始出现问题...

行列式应该是N+1。

计算U对角线的乘积可能是误差的传播：

from scipy.sparse import diags
from scipy.sparse.linalg import splu
from sksparse.cholmod import cholesky
from math import exp

n=int(5e6)
K = diags([-1.],-1,shape=(n,n)) + diags([2.],shape=(n,n)) + diags([-1.],1,shape=(n,n))
lu = splu(K.tocsc())
diagL = lu.L.diagonal()
diagU = lu.U.diagonal()
det=diagL.prod()*diagU.prod()
print(det)

factor = cholesky(K.tocsc())
ld = factor.logdet()
print(exp(ld))

输出：

4999993.625461911

4999993.625461119

即使 U 精确到 10-13 位，这也是可以预期的：

n=int(5e6)
print(n*diags([1-0.00000000000025],0,shape=(n,n)).diagonal().prod())

4999993.749444371