在 Julia 中计算稀疏矩阵的对数行列式
Computing the logdeterminant of a Sparse Matrix in Julia
我对计算大型、稀疏、复杂(浮点)矩阵的行列式的对数很感兴趣。我的第一个想法是使用 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
将始终打印出 "NOT correct" 选项。此外,
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 安装是否有问题(尽管我已经在另一台电脑)。如果您能给我任何建议,那就太好了!
原因是稀疏 LU 也有一个比例因子,可以用 LUFs[:Rs]
提取(在 Julia 0.6 中和 LUFs.Rs
在 Julia 0.7- 中)。因此计算变为
julia> det(LUFs[:U])/prod(LUFs[:Rs])
0.4576579970452131 - 0.07585833005688908im
julia> det(A)
0.4576579970452133 - 0.07585833005688908im
对于稀疏情况,我们可能也应该有 logabsdet
。但是,您的矩阵是正定的吗?这样,您就可以计算 Cholesky 分解的 logdet
。
我对计算大型、稀疏、复杂(浮点)矩阵的行列式的对数很感兴趣。我的第一个想法是使用 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
将始终打印出 "NOT correct" 选项。此外,
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 安装是否有问题(尽管我已经在另一台电脑)。如果您能给我任何建议,那就太好了!
原因是稀疏 LU 也有一个比例因子,可以用 LUFs[:Rs]
提取(在 Julia 0.6 中和 LUFs.Rs
在 Julia 0.7- 中)。因此计算变为
julia> det(LUFs[:U])/prod(LUFs[:Rs])
0.4576579970452131 - 0.07585833005688908im
julia> det(A)
0.4576579970452133 - 0.07585833005688908im
对于稀疏情况,我们可能也应该有 logabsdet
。但是,您的矩阵是正定的吗?这样,您就可以计算 Cholesky 分解的 logdet
。