25

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

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

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

4

4 回答 4

7

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

  • 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 等。等人，2012]。主要技巧在于我们可以将 log(det(Σ)) 写为 trace(log(Σ))，其中 log(Σ) 是矩阵对数。

于 2014-02-19T06:11:32.040 回答
6

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

于 2013-10-27T10:23:00.567 回答
5

您可以用于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 = 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
于 2020-04-01T23:41:17.590 回答
1

使用 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

于 2021-04-27T18:22:04.127 回答