[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