[Numpy-svn] r8384 - in trunk/numpy/linalg: . tests
numpy-svn at scipy.org
numpy-svn at scipy.org
Tue May 4 23:34:47 EDT 2010
Author: charris
Date: 2010-05-04 22:34:47 -0500 (Tue, 04 May 2010)
New Revision: 8384
Modified:
trunk/numpy/linalg/__init__.py
trunk/numpy/linalg/linalg.py
trunk/numpy/linalg/tests/test_linalg.py
Log:
ENH: Add slogdet to the linalg module. The patch is from njs with
slogdet substituted for sign_log_det. Closes ticket #1402.
Modified: trunk/numpy/linalg/__init__.py
===================================================================
--- trunk/numpy/linalg/__init__.py 2010-05-05 02:50:03 UTC (rev 8383)
+++ trunk/numpy/linalg/__init__.py 2010-05-05 03:34:47 UTC (rev 8384)
@@ -9,6 +9,7 @@
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
Modified: trunk/numpy/linalg/linalg.py
===================================================================
--- trunk/numpy/linalg/linalg.py 2010-05-05 02:50:03 UTC (rev 8383)
+++ trunk/numpy/linalg/linalg.py 2010-05-05 03:34:47 UTC (rev 8384)
@@ -10,15 +10,15 @@
"""
__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,10 +1533,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
----------
a : array_like, shape (M, M)
@@ -1544,9 +1549,16 @@
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
-----
The determinant is computed via LU factorization using the LAPACK
@@ -1557,9 +1569,23 @@
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,11 +1603,51 @@
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
def lstsq(a, b, rcond=-1):
Modified: trunk/numpy/linalg/tests/test_linalg.py
===================================================================
--- trunk/numpy/linalg/tests/test_linalg.py 2010-05-05 02:50:03 UTC (rev 8383)
+++ trunk/numpy/linalg/tests/test_linalg.py 2010-05-05 03:34:47 UTC (rev 8384)
@@ -127,13 +127,32 @@
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):
u, s, vt = linalg.svd(a, 0)
More information about the Numpy-svn
mailing list