Log Determinant?

[John Myles White]

Do we have a log determinant function? I'm in the process of implementing a quick one out of necessity, but would prefer to work with something built into Julia.

-- John

[Tim Holy]

Is there something better than Cholesky/LU followed by product of log of
diagonals? If not, that would be a few-liner and I think should probably be in
Base, it's just so gosh-durnit useful.

--Tim

[John Myles White]

After some web-searching I couldn't find anything better. Right now I'm using the one-liner:

logdet(A::Matrix) = sum(log(diag(chol(A))))

I'll add that to Base now and let others replace it if there's something better available.

-- John

[Tim Holy]

On Tuesday, February 26, 2013 09:32:47 PM John Myles White wrote:
> After some web-searching I couldn't find anything better. Right now I'm
> using the one-liner:
>
> logdet(A::Matrix) = sum(log(diag(chol(A))))
>
> I'll add that to Base now and let others replace it if there's something
> better available.

You can only use chol if it's symmetric positive definite. The vast majority of
practical uses of this are, but you need an LU-based fallback in cases where
that's not true. (The matrix can even have negative eigenvalues as long as
there are an even number of them, and logdet would still be well-defined.)

--Tim

[John Myles White]

Should I add an explicit test for symmetry or just always use an LU-decomposition?

-- John

[Tim Holy]

On Tuesday, February 26, 2013 09:54:32 PM John Myles White wrote:
> Should I add an explicit test for symmetry or just always use an
> LU-decomposition?

Or perhaps a two-argument version, logdet(A::Matrix, isposdef::Bool). Perhaps
if the second argument is omitted, default to the LU version.

[Simon Byrne]

Perhaps it should be logabsdet, since (I assume) you'll be discarding the sign from the LU? 

> logdet(A::Matrix) = sum(log(diag(chol(A)))) 

Also, I think you mean

logdet(A::Matrix) = 2*sum(log(diag(chol(A)))) 

(I've made the same mistake before: it took me a week to find the cause).

You could also make LU the default for Matrix, and a different method for the factorization object?

[John Myles White]

On Feb 27, 2013, at 7:13 AM, Simon Byrne <simon...@gmail.com> wrote:

Perhaps it should be logabsdet, since (I assume) you'll be discarding the sign from the LU? 

> logdet(A::Matrix) = sum(log(diag(chol(A)))) 

For my purposes discarding the sign is fine, but I'm not sure for others.

Also, I think you mean

logdet(A::Matrix) = 2*sum(log(diag(chol(A)))) 

(I've made the same mistake before: it took me a week to find the cause).

Good catch. I fixed this before copying the code to Base.

You could also make LU the default for Matrix, and a different method for the factorization object?

This seems reasonable as well. I think Tim's suggestion seemed good because you typically want to do this in a case where you don't have a factorization, but do know that the matrix admits the chol() call.

I feel like the more linear algebraic folks here should make the final call on API.

 -- John

[Steven G. Johnson]

On Wednesday, February 27, 2013 7:13:46 AM UTC-5, Simon Byrne wrote:

Perhaps it should be logabsdet, since (I assume) you'll be discarding the sign from the LU? 

In the applications where I run into logdet (mostly in quantum field theory), the matrix is not necessarily Hermitian positive-definite (so you have to use LU not Cholesky), but for physical (and algebraic) reasons the determinant has to be positive.

If the determinant were negative, I wouldn't want Julia to silently discard the sign, I would want it to throw a DomainError() exception, since it would indicate an error on my part.

--SGJ

[Douglas Bates]

The convention in R is to return the logarithm of the absolute value of the determinant and separately return the sign, which is easy to calculate from the permutation vector.