[Numpy-discussion] Best interface for computing the logarithm of the determinant?
Nathaniel Smith
njs@pobox....
Fri Feb 19 02:08:50 CST 2010
Thanks for your comments, all. Since it occurs to me that this is a
general need, not just for sparse matrices, and it would be very
annoying to settle on one API for scikits.sparse and then have another
show up in one of the main packages later, I've just submitted a patch
for option (1) to numpy:
http://projects.scipy.org/numpy/ticket/1402
And we'll see what happens with that :-)
On Tue, Feb 16, 2010 at 12:49 PM, Nathaniel Smith <njs@pobox.com> 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
>
> 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
>
> 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)
>
> These are all kind of ugly looking, unfortunately, but that seems
> unavoidable, unless someone has a clever idea.
>
> Any preferences?
>
> -- Nathaniel
>
More information about the NumPy-Discussion
mailing list