[SciPy-dev] Best interface for computing the logarithm of the determinant?
Anne Archibald
peridot.faceted@gmail....
Tue Feb 16 17:14:02 CST 2010
On 16 February 2010 16:19, Dag Sverre Seljebotn
<dagss@student.matnat.uio.no> wrote:
> 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)
I kind of like this one, but it could definitely become a pain. For
what it's worth, for complex numbers it means you get a polar
representation of the complex number (complex number of magnitude 1
and real log of the magnitude).
>> 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).
I don't really like this very much, since it can be really startling
to suddenly have complex numbers appear out of nowhere. But at least
it's consistent.
>>
>> 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.
This lends itself well to a multiply-and-rescale-when-necessary
approach, and avoids transcendental functions entirely; if the 10
offends, you could use 2 or 16, I suppose. But this kind of
count-every-flop approach is not really appropriate for scipy/numpy.
>> 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.
Since dets are (potentially) expensive operations, no, I don't like
this either. Of course, one might want a sgndet function in addition
to the usual API.
Anne
More information about the SciPy-Dev
mailing list