[SciPy-dev] Best interface for computing the logarithm of the determinant?
josef.pktd@gmai...
josef.pktd@gmai...
Tue Feb 16 15:29:12 CST 2010
On Tue, Feb 16, 2010 at 4:19 PM, 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)
>>
>> 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'm also in favor of this, but would prefer if promotion to complex
only occurs if the result requires it.
np.real_if_close might be safer than .real
Josef
>
>>
>> 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
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
More information about the SciPy-Dev
mailing list