下面是我在回答here时提供的一些参考，我认为它们解决了您试图解决的实际问题：

• notes，用于幕府图书馆中的实现
• 埃尔兰德·奥恩，丹尼尔·辛普森：* 高维高斯分布中的参数估计 *，特别是第2.1节（arxiv:1105.5256）
• 伊尔斯·C.F.易普森、迪安·J·李：* 行列式近似 *（arxiv:1105.0437）
• 阿诺德·鲁斯肯：* 大型稀疏对称正定矩阵行列式的逼近 *（arxiv:hep-lat/0008007）

引用幕府笔记：
计算似然表达式中的对数行列式项的常用技术依赖于矩阵的Cholesky因式分解，即Σ=LLT，（L是下三角Cholesky因子），然后使用该因子的对角项计算log（检测（Σ））=2∑ni= 1 log然而，对于稀疏矩阵，如协方差矩阵通常那样，Cholesky因子经常遭受填充现象-它们本身并不那么稀疏。因此，对于大的维数，由于需要大量的存储器来存储因子的所有这些不相关的非对角系数，所以该技术变得不可行。虽然已经开发了排序技术来预先置换行和列以减少填充，例如近似最小度（AMD）重新排序，但是这些技术很大程度上依赖于稀疏模式，因此不能保证给予更好的结果。
最近的研究表明，使用复分析、数值线性代数和贪婪图着色等技术，我们可以将对数行列式逼近到任意精度[Aune et. al.，2012]。主要技巧在于我们可以将log（det（Σ））写成trace（log（Σ）），其中log（Σ）是矩阵对数。

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

可以使用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 = (-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

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