[SciPy-dev] Best interface for computing the logarithm of the determinant?

Dag Sverre Seljebotn dagss at student.matnat.uio.no
Tue Feb 16 16:19:58 EST 2010


Nathaniel Smith wrote:
> So when you have a matrix whose determinant you want, it's often wise
> to compute the logarithm of the determinant instead of the determinant
> itself, because determinants involve lots and lots of multiplications
> and the result might otherwise underflow/overflow. Therefore, in
> scikits.sparse, I'd like to provide an API for doing this (and this is
> well-supported by the underlying libraries).
> 
> But the problem is that for a general matrix, the determinant may be
> zero or negative. Obviously we can deal with this, but what's the best
> API? I'd like to use one consistently across the different
> factorizations in scikits.sparse, and perhaps eventually in numpy as
> well.
> 
> Some options:
> 
> 1) Split off the sign into a separate return value ('sign' may be 1, -1, 0):
> sign, value = logdet(A)
> actual_determinant = sign * exp(value)
> 
> 2) Allow complex/infinite return values, even when A is a real matrix:
> logdet(eye(3)) == pi*1j
> logdet(zeros((3, 3))) == -Inf

I'm +1 for this one. It is easy to interpret, and if one knows that the 
result is positive (like after a successful CHOLMOD LL^T) is is easy 
enough to take the real part by appending .real.

Also it makes more sense for code that may be dealing with both complex 
and real matrices compared to 1).

> 
> 3) "Scientific notation" (This is what UMFPACK's API does): return a
> mantissa and base-10 exponent:
> mantissa, exponent = logdet(A)
> actual_determinant = mantissa * 10 ** exponent

This is my least preferred option, because "10" has nothing to do with 
floats.

> 
> 4) Have separate functions for computing the sign, and the log of the
> absolute value (This is what GSL does, though it seems pointlessly
> inefficient):
> sign = sgndet(A)
> value = logdet(A)
> actual_determinant = sign * exp(value)

Well, it doesn't have to be inefficient, you use cache the return value 
internally on id(A) and possibly weak reference callbacks to remove the 
cached values...but pretty ugly, yeah.

-- 
Dag Sverre



More information about the SciPy-Dev mailing list