CODEWITHSUNDEEP
juliasparse-matrixdeterminantsmatrix-factorization
MDM
MDM

Reputation: 65

Computing the logdeterminant of a Sparse Matrix in Julia

I'm interested in computing the log of the determinant of a large, sparse, complex (floating point) matrix. My first thought is to use an LU decomposition, i.e.:

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

which will always print out the "NOT correct" option. Furthermore, 

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

always comes out as false.

If instead I try to directly use

logdet(LUFs)

or

logdet(sparse(A))

I get an error:

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

I'm not sure if it's something wrong in the way I've coded it (I'm a beginner transitioning from matlab), or if there's something wrong with my Julia install (though I have replicated these results on another computer). Any pointers you could give me would be great!

Upvotes: 2

Views: 452

Answers (1)

Andreas Noack
Andreas Noack

Reputation: 1380

The reason is that the sparse LU has a scaling factor as well which can be extracted with LUFs[:Rs] (in Julia 0.6 and LUFs.Rs in Julia 0.7-). Hence the computation becomes

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

julia> det(A)
0.4576579970452133 - 0.07585833005688908im

We should probably have logabsdet for the sparse case as well. However, is your matrix positive definite by any chance? In so, you'd be able to compute the logdet of the Cholesky factorization.

Upvotes: 2

Related Questions