在 Julia 中计算稀疏矩阵的对数行列式答案

【问题标题】：Computing the logdeterminant of a Sparse Matrix in Julia在 Julia 中计算稀疏矩阵的对数行列式
【发布时间】：2018-06-03 09:45:32
【问题描述】：

我有兴趣计算大型、稀疏、复杂（浮点）矩阵的行列式的对数。我的第一个想法是使用 LU 分解，即：

srand(123) 
A=complex.(rand(3,3), rand(3,3))

LUF=lufact(A)
LUFs=lufact(sparse(A))

if round(det(LUFs[:L])*det(LUFs[:U])-det(A[LUFs[:p], LUFs[:q]]), 5)==0
    println("Sparse LU determinant is correct\n")
else
    println("Sparse LU determinant is NOT correct\n")
end

这将始终打印出“不正确”选项。此外，

round.(LUFs[:L]*LUFs[:U], 5)==round.(A[LUFs[:p], LUFs[:q]], 5)

总是假的。

如果我尝试直接使用

logdet(LUFs)

或

logdet(sparse(A))

我收到一个错误：

LoadError: MethodError: no method matching
logabsdet(::Base.SparseArrays.UMFPACK.UmfpackLU{Complex{Float64},Int64})
Closest candidates are:
logabsdet(!Matched::Base.LinAlg.UnitUpperTriangular{T,S} where S<:
(AbstractArray{T,2} where T)) where T at linalg/triangular.jl:2184
logabsdet(!Matched::Base.LinAlg.UnitLowerTriangular{T,S} where S<:
(AbstractArray{T,2} where T)) where T at linalg/triangular.jl:2185
logabsdet(!Matched::Union{LowerTriangular{T,S} where S<:
(AbstractArray{T,2} where T), UpperTriangular{T,S} where S<:
(AbstractArray{T,2} where T)}) where T at linalg/triangular.jl:2189
...
while loading...

我不确定我的编码方式是否有问题（我是从 matlab 过渡的初学者），或者我的 Julia 安装是否有问题（尽管我已经在另一台计算机上复制了这些结果）。你能给我的任何指点都会很棒！

【问题讨论】：

距离有多近？这可能只是浮点舍入错误。
您几乎不应该使用== 来比较浮点计算的结果，尤其是通过多线程 BLAS 的结果。这适用于任何语言的编程。由于这些数字的工作方式，会有一些舍入误差。
这是一个非常好的观点，但它偏离了一个数量级。对于上面的 rng 种子，我得到 det(LUFs[:L])*det(LUFs[:U])=0.0216… - 0.00358…im 和 det(A)=0.458… - 0.0759…im。
（针对这一点进行了编辑）
速度有多重要？使用 SVD 分解并执行类似log(complex(det(U))) + sum(log.(d)) + logdet(complex(det(T))) 的操作。

标签： julia sparse-matrix determinants matrix-factorization

【解决方案1】：

原因是稀疏 LU 也有一个缩放因子，可以用 LUFs[:Rs] 提取（在 Julia 0.6 中和在 Julia 0.7- 中的 LUFs.Rs）。因此计算变成了

julia> det(LUFs[:U])/prod(LUFs[:Rs])
0.4576579970452131 - 0.07585833005688908im

julia> det(A)
0.4576579970452133 - 0.07585833005688908im

对于稀疏的情况，我们可能也应该有logabsdet。但是，您的矩阵是否是正定的？这样，您就可以计算 Cholesky 分解的logdet。

【讨论】：

不幸的是我的矩阵不是正定的，但是添加比例因子可以解决我的问题。谢谢！