Provide logabsdet(::Matrix)?

[bsxfan]

The current logdet implementation in linalg/dense.jl,

logdet(A::Matrix) = 2.0 * sum(log(diag(cholfact(A)[:U])))

would probably be cleaner if implemented as:

logdet(A::Matrix) = logdet(cholfact(A))

[JeffBezanson]

There is a bug in cholfact:

julia> logdet(rand(3,3))
ERROR: PosDefException not defined
 in cholfact! at linalg/factorization.jl:13
 in logdet at linalg/dense.jl:422

Also, is there or can there be some interface for getting log(abs(det)) and the sign? Or det of non-positive-definite matrices?

[dmbates]

@JeffBezanson A det method for the LU type already exists. So logdet for a square matrix should be defined similarly although it would need to return log(abs(det)) and the sign separately.

[bsxfan]

det() and LU are already equipped to return a complex result for a complex matrix. So maybe logdet(::LU) could always return a complex value too, which would accomodate the case of both negative and complex determinants.

The inconsistency is that at present logdet(::Matrix) returns a real result. Maybe one should rather define complex_logdet() which always returns a complex result.

[StefanKarpinski]

Having logdet always return a complex result occurred to me too and feels more appealing than returning both a value and a sign, bu it does seem inconsistent with our other functions, which return NaN when applied to a real value in such a way that the result would be non-real.

[ViralBShah]

Cc: @dmbates @andreasnoack

[JeffBezanson]

We now throw domain errors instead of returning NaN.

[StefanKarpinski]

Oh, right. But same idea.

[dmbates]

The thinking among the R developers was that returning the logarithm of the absolute value of the determinant along with its sign was done to provide a more stable calculation and to avoid overflow/underflow. However I will admit that almost all the time I use logdet it is for a positive definite matrix, usually some form of a variance-covariance matrix.

[bsxfan]

While we're on logdet, it might be nice to supply Woodbury with a logdet too. Woodbury is e.g. useful in probability theory when doing factor analysis, in which case the determinants also end up being positive. Since the woodbury matrices can be huge, the determinant may well under/overflow, but logdet won´t.

[stevengj]

Restricting logdet to Hermitian positive-definite matrices seems wrong to me.

For example, in some problems arising in integral-equation methods for quantum field theory, we need the logdet of a real matrix A that is nonsymmetric but its symmetric part is positive definite, i.e. A + A' is symmetric positive-definite. Such a matrix always has a positive determinant and so its logdet is real.

Also, currently the logdet of a Hermitian positive-definite matrix returns a complex result (with zero imaginary part), which doesn't seem right.

[dmbates]

I agree with @stevengj that returning a complex result from logdet of a Hermitian positive-definite matrix doesn't seem right.

[StefanKarpinski]

So returning a complex number from logdet doesn't seem right and restricting logdet to positive-definite matrices doesn't seem right either. I can see only two other options:

1. Return both a real value and a sign flag from logdet.
2. Restrict logdet to some larger set of matrix inputs – possibly any matrix for which the determinant is positive?

Are there other options I'm missing? If not, which of these is the better choice?

[andreasnoack]

It could work like eig and return a complex result when the determinant is negative and real otherwise.

[StefanKarpinski]

I'm not a fan of eig behaving that way, honestly.

[andreasnoack]

Then I'd prefer DomainError over tuple output.

[timholy]

Another option is to have two functions, of which one could be called logabsdet for example.

[stevengj]

@StefanKarpinski, I would tend to (2) restrict logdet to matrices with positive real determinants, in which case it always returns a real result, and throw a DomainError otherwise. To begin with, there are three main cases:

• positive-definite (either real-symmetric or Hermitian): use Cholesky factorization
• real but unsymmetric (but with positive determinant): use LU
• real-symmetric/Hermitian but indefinite (with positive determinant): use LDL* factorization

In theory, I suppose that you could also have an unsymmetrical complex matrix that nevertheless has a positive real determinant, but I'm not sure that this case arises enough that it is worth the trouble to figure out how to compute logdet efficiently for it.

[ViralBShah]

How about having logdet as proposed by @stevengj and logdetabs, a slight modification of the name proposed by @timholy, to make the name more discoverable?

[dmbates]

@ViralBShah, @timholy Tim proposed logabsdet not logdetabs. I think the logabsdet name more clearly reflects the operation being performed.

[pao]

I think Viral was acknowledging that, but logdetabs is tab-completable from logdet.

[StefanKarpinski]

Another option would be a keyword argument indicating whether to take the absolute value or not.

[ViralBShah]

Of course - let's go with a keyword argument.

[stevengj]

Could someone enlighten me on the advantage of logdet(A, abs=true) over logdetabs(A)?

[ViralBShah]

I was proposing logdet(A; abs=true). The real question is whether the two routines are similar enough that they should share a name.

[timholy]

I think @stevengj 's point was that logdetabs is shorter to type than the keyword/optional input version.

But if we do give it a separate name, I think logabsdet is really much better than logdetabs, despite the tab-completion issue. logdetabs(A) = log(det(abs(A))), but logabsdet(A) = log(abs(det(A))). These are very different, and it would be dangerous to confuse them.

[StefanKarpinski]

Yeah, I'm not sold on using a keyword, just an idea. Calling it logdetabs seems outright obtuse.

[mlubin]

It sounds like logabsdet and logdet return different things. Is it good for keyword arguments to change the function output type?

[StefanKarpinski]

Yeah, good point. I was imagining the keyword as only affecting whether a negative determinant causes an error or a sign flip (I realize it may not be implemented that way).

[ViralBShah]

Ok. The arguments for logabsdet are quite persuasive.

[bsxfan]

See my comment in #16, where I give a tyre-kicking implementation of a generalized, but type-stable logdet for the LU factorization. There is a logdet2 which returns log(abs(det)) and sign(det) as well as a logdet which returns complex log for complex arguments and real log for real arguments. The latter gives Domainerror when the det<0, but works for non-positive definite matrices with positive determinants, e.g. diagm(-ones(2)).

[garrison]

FWIW, in numpy the logabsdet() function is called slogdet()