2

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!

MDM
  • 65
  • 1
  • 6
  • 1
    How close is it? It might just be floatingpoint rounding error. – Oscar Smith Dec 20 '17 at 19:11
  • 1
    You should almost never use `==` to compare the results of a floating point calculation, especially ones that go through a multithreaded BLAS. This is true for programming in any language. There will be some rounding error just by virtue of how these numbers work. – Chris Rackauckas Dec 20 '17 at 19:47
  • That's a really good point, but it's off by an order of magnitude. For the rng seed above, I get det(LUFs[:L])*det(LUFs[:U])=0.0216… - 0.00358…im and det(A)=0.458… - 0.0759…im. – MDM Dec 21 '17 at 00:32
  • (edited to address this point) – MDM Dec 21 '17 at 00:39
  • How much does speed matter? Use SVD decomposition and do something like `log(complex(det(U))) + sum(log.(d)) + logdet(complex(det(T)))`. – carstenbauer Dec 21 '17 at 09:01
  • speed is probably somewhat flexible, but trying `svdfact(sparse(A))`, I get a similar error message to that for `logdet(sparse(A))`. – MDM Dec 21 '17 at 19:26

1 Answers1

2

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.

Andreas Noack
  • 1,370
  • 6
  • 11
  • Unfortunately my matrix isn't positive definite, but adding in the scaling factor fixes my issue. Thank you! – MDM Dec 21 '17 at 19:13