[WIP] add logabsdet(::UmfpackLU) method

[ElOceanografo]

Will fix JuliaLang/LinearAlgebra.jl#788. This new method is a direct translation of the dense logabsdet(::LU) method in LinearAlgebra:

julia/stdlib/LinearAlgebra/src/lu.jl

Line 457 in 3365b4e

function logabsdet(F::LU{T}) where T # return log(abs(det)) and sign(det)

It produces the correct answers, but is relatively slow and needs work on optimization. The main bottleneck is extracting the factors via umf_extract, so directly indexing into the UmfpackLU struct should speed it up...though I'm not sure how to do that.

[ViralBShah]

@ElOceanografo Maybe address the comments and get this merged?

[ElOceanografo]

Sorry, lost track of this! Comments were straightforward, but this is still very slow. Benchmarking on my computer produces the following timings:

using LinearAlgebra, SparseArrays, BenchmarkTools
n = 2000
pnz = 0.01
A = sprandn(n, n, pnz); A = A'A + I
Adense = Matrix(A)
F = lu(A)
Fdense = lu(Adense)

@btime logabsdet($A) # 549.719 ms (91 allocations: 278.57 MiB)
@btime logabsdet($F) # 79.753 ms (25 allocations: 86.26 MiB)
@btime logabsdet($Adense) # 46.192 ms (3 allocations: 30.53 MiB)
@btime logabsdet($Fdense) # 20.099 μs (0 allocations: 0 bytes)

Profiling shows most of the time in the first benchmark is spent here, while most of the time in the second benchmark is spent extracting the fields from the LU factorization in this line. I know next to nothing about the internals of SuiteSparse, maybe someone who does can chime in? Using the sparse factorization should in principle be faster than converting to a dense matrix and falling back on the default method, right?

Edit: added timing for logabsdet(Fdense)

[ViralBShah]

The problem with testing on random sparse matrices is that the LU factorization gets a lot of fill-in, making it almost dense. In this case, the filled up factors end up having 50% nonzeros. In reality, most sparse problems coming from a real application have more structure and significantly less fill-in.

[ViralBShah]

UmfpackLU is an opaque object, but this ends up spending a bunch of time converting the opaque object, when all you really want are the diagonals. I still think this should be ok to merge. Try some experiments with the identity matrix?

At scale, the U[i,i] indexing could end up being slow, but it is easy to pick that entry from the sparse data structure directly, much faster. Worth doing some experiments with a million size matrix.

[ViralBShah]

Probably also rebase to master.

[ElOceanografo]

Running the same benchmark with A = spdiagm(0 => ones(n)) is still slower than the dense version (189.999 μs (81 allocations: 1.11 MiB)). If you know where the UmfpackLU object's layout is documented, I can try to figure out how to pick the diagonal values out directly. But if you think this is worth merging as-is I can go ahead and rebase it.

[ViralBShah]

That will have to be in the UmfpackLU documentation or maybe even the Umfpack source code. I had trouble compiling this PR (following the gh command line) instructions and did end up rebasing it to master. Trying to get it working so I can try a few examples myself.

[ElOceanografo]

I did some digging in the Umfpack documentation, and it should be possible to extract just the U matrix and scaling factors Rs by passing null values in for the other arguments (see section 13, top of page 94, in the pdf here: https://fossies.org/linux/SuiteSparse/UMFPACK/Doc/UMFPACK_UserGuide.pdf). That isn't currently possible through the SuiteSparse.UMFPACK wrappers, which is probably worth a separate issue.

[ViralBShah]

Something has gone wrong with this PR.

[ElOceanografo]

I tried to rebase it to master last night and may have messed it up. I can try to fix it, or if it's easier I could close it and open a new one with the same changes...

[ViralBShah]

If its easier to do a new PR, just go for it.

[ElOceanografo]

Ok, closing this, will open a new one in a second.