ENH: Add slogdet to the linalg module. The patch is from njs with

slogdet substituted for sign_log_det. Closes ticket #1402.

3 files changed, 97 insertions, 11 deletions

diff --git a/numpy/linalg/__init__.py b/numpy/linalg/__init__.py
index 58b012f00..a74a31950 100644
--- a/numpy/linalg/__init__.py
+++ b/numpy/linalg/__init__.py
@@ -9,6 +9,7 @@ norm            Vector or matrix norm
 inv             Inverse of a square matrix
 solve           Solve a linear system of equations
 det             Determinant of a square matrix
+slogdet         Logarithm of the determinant of a square matrix
 lstsq           Solve linear least-squares problem
 pinv            Pseudo-inverse (Moore-Penrose) calculated using a singular
                 value decomposition
diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py
index 005738c1c..a0d191b4b 100644
--- a/numpy/linalg/linalg.py
+++ b/numpy/linalg/linalg.py
@@ -10,15 +10,15 @@ zgetrf, dpotrf, zpotrf, dgeqrf, zgeqrf, zungqr, dorgqr.
 """
 
 __all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv',
-           'cholesky', 'eigvals', 'eigvalsh', 'pinv', 'det', 'svd',
-           'eig', 'eigh','lstsq', 'norm', 'qr', 'cond', 'matrix_rank',
+           'cholesky', 'eigvals', 'eigvalsh', 'pinv', 'slogdet', 'det',
+           'svd', 'eig', 'eigh','lstsq', 'norm', 'qr', 'cond', 'matrix_rank',
            'LinAlgError']
 
 from numpy.core import array, asarray, zeros, empty, transpose, \
         intc, single, double, csingle, cdouble, inexact, complexfloating, \
         newaxis, ravel, all, Inf, dot, add, multiply, identity, sqrt, \
         maximum, flatnonzero, diagonal, arange, fastCopyAndTranspose, sum, \
-        isfinite, size, finfo
+        isfinite, size, finfo, absolute, log, exp
 from numpy.lib import triu
 from numpy.linalg import lapack_lite
 from numpy.matrixlib.defmatrix import matrix_power
@@ -1533,9 +1533,14 @@ def pinv(a, rcond=1e-15 ):
 
 # Determinant
 
-def det(a):
+def slogdet(a):
     """
-    Compute the determinant of an array.
+    Compute the sign and (natural) logarithm of the determinant of an array.
+
+    If an array has a very small or very large determinant, than a call to
+    `det` may overflow or underflow. This routine is more robust against such
+    issues, because it computes the logarithm of the determinant rather than
+    the determinant itself.
 
     Parameters
     ----------
@@ -1544,8 +1549,15 @@ def det(a):
 
     Returns
     -------
-    det : ndarray
-        Determinant of `a`.
+    sign : float or complex
+        A number representing the sign of the determinant. For a real matrix,
+        this is 1, 0, or -1. For a complex matrix, this is a complex number
+        with absolute value 1 (i.e., it is on the unit circle), or else 0.
+    logdet : float
+        The natural log of the absolute value of the determinant.
+
+    If the determinant is zero, then `sign` will be 0 and `logdet` will be
+    -Inf. In all cases, the determinant is equal to `sign * np.exp(logdet)`.
 
     Notes
     -----
@@ -1557,9 +1569,23 @@ def det(a):
     The determinant of a 2-D array [[a, b], [c, d]] is ad - bc:
 
     >>> a = np.array([[1, 2], [3, 4]])
-    >>> np.linalg.det(a)
+    >>> (sign, logdet) = np.linalg.slogdet(a)
+    >>> (sign, logdet)
+    (-1, 0.69314718055994529)
+    >>> sign * np.exp(logdet)
     -2.0
 
+    This routine succeeds where ordinary `det` does not:
+
+    >>> np.linalg.det(np.eye(500) * 0.1)
+    0.0
+    >>> np.linalg.slogdet(np.eye(500) * 0.1)
+    (1, -1151.2925464970228)
+
+    See Also
+    --------
+    det
+
     """
     a = asarray(a)
     _assertRank2(a)
@@ -1577,10 +1603,50 @@ def det(a):
     if (info < 0):
         raise TypeError, "Illegal input to Fortran routine"
     elif (info > 0):
-        return 0.0
-    sign = add.reduce(pivots != arange(1, n+1)) % 2
-    return (1.-2.*sign)*multiply.reduce(diagonal(a), axis=-1)
+        return (t(0.0), _realType(t)(-Inf))
+    sign = 1. - 2. * (add.reduce(pivots != arange(1, n + 1)) % 2)
+    d = diagonal(a)
+    absd = absolute(d)
+    sign *= multiply.reduce(d / absd)
+    log(absd, absd)
+    logdet = add.reduce(absd, axis=-1)
+    return sign, logdet
+
+def det(a):
+    """
+    Compute the determinant of an array.
 
+    Parameters
+    ----------
+    a : array_like, shape (M, M)
+        Input array.
+
+    Returns
+    -------
+    det : ndarray
+        Determinant of `a`.
+
+    Notes
+    -----
+    The determinant is computed via LU factorization using the LAPACK
+    routine z/dgetrf.
+
+    Examples
+    --------
+    The determinant of a 2-D array [[a, b], [c, d]] is ad - bc:
+
+    >>> a = np.array([[1, 2], [3, 4]])
+    >>> np.linalg.det(a)
+    -2.0
+
+    See Also
+    --------
+    slogdet : Another way to representing the determinant, more suitable
+      for large matrices where underflow/overflow may occur.
+
+    """
+    sign, logdet = slogdet(a)
+    return sign * exp(logdet)
 
 # Linear Least Squares
 
diff --git a/numpy/linalg/tests/test_linalg.py b/numpy/linalg/tests/test_linalg.py
index b81561254..3a9584dd6 100644
--- a/numpy/linalg/tests/test_linalg.py
+++ b/numpy/linalg/tests/test_linalg.py
@@ -127,12 +127,31 @@ class TestPinv(LinalgTestCase, TestCase):
 class TestDet(LinalgTestCase, TestCase):
     def do(self, a, b):
         d = linalg.det(a)
+        (s, ld) = linalg.slogdet(a)
         if asarray(a).dtype.type in (single, double):
             ad = asarray(a).astype(double)
         else:
             ad = asarray(a).astype(cdouble)
         ev = linalg.eigvals(ad)
         assert_almost_equal(d, multiply.reduce(ev))
+        assert_almost_equal(s * np.exp(ld), multiply.reduce(ev))
+        if s != 0:
+            assert_almost_equal(np.abs(s), 1)
+        else:
+            assert_equal(ld, -inf)
+
+    def test_zero(self):
+        assert_equal(linalg.det([[0.0]]), 0.0)
+        assert_equal(type(linalg.det([[0.0]])), double)
+        assert_equal(linalg.det([[0.0j]]), 0.0)
+        assert_equal(type(linalg.det([[0.0j]])), cdouble)
+
+        assert_equal(linalg.slogdet([[0.0]]), (0.0, -inf))
+        assert_equal(type(linalg.slogdet([[0.0]])[0]), double)
+        assert_equal(type(linalg.slogdet([[0.0]])[1]), double)
+        assert_equal(linalg.slogdet([[0.0j]]), (0.0j, -inf))
+        assert_equal(type(linalg.slogdet([[0.0j]])[0]), cdouble)
+        assert_equal(type(linalg.slogdet([[0.0j]])[1]), double)
 
 class TestLstsq(LinalgTestCase, TestCase):
     def do(self, a, b):