Enable doctest ELLIPSIS in a few modules.

[dunshire.git] / dunshire / matrices.py

diff --git a/dunshire/matrices.py b/dunshire/matrices.py
index 6ed4f851487235dba117ef1521151bd432db9076..94f841d39f888392624077afd12bcdcb7e0c3c19 100644 (file)
--- a/dunshire/matrices.py
+++ b/dunshire/matrices.py
@@ -3,9 +3,10 @@ Utility functions for working with CVXOPT matrices (instances of the
 class:`cvxopt.base.matrix` class).
 """
 
+from copy import copy
 from math import sqrt
 from cvxopt import matrix
-from cvxopt.lapack import gees, syevr
+from cvxopt.lapack import gees, gesdd, syevr
 
 from . import options
 
@@ -140,7 +141,10 @@ def eigenvalues(symmat):
 
     domain_dim = symmat.size[0]
     eigs = matrix(0, (domain_dim, 1), tc='d')
-    syevr(symmat, eigs)
+
+    # Create a copy of ``symmat`` here because ``syevr`` clobbers it.
+    dummy = copy(symmat)
+    syevr(dummy, eigs)
     return list(eigs)
 
 
@@ -209,7 +213,7 @@ def eigenvalues_re(anymat):
     return [eig.real for eig in eigs]
 
 
-def identity(domain_dim):
+def identity(domain_dim, typecode='i'):
     """
     Create an identity matrix of the given dimensions.
 
@@ -219,6 +223,10 @@ def identity(domain_dim):
     domain_dim : int
         The dimension of the vector space on which the identity will act.
 
+    typecode : {'i', 'd', 'z'}, optional
+        The type code for the returned matrix, defaults to 'i' for integers.
+        Can also be 'd' for real double, or 'z' for complex double.
+
     Returns
     -------
 
@@ -247,7 +255,7 @@ def identity(domain_dim):
     entries = [int(i == j)
                for i in range(domain_dim)
                for j in range(domain_dim)]
-    return matrix(entries, (domain_dim, domain_dim))
+    return matrix(entries, (domain_dim, domain_dim), tc=typecode)
 
 
 def inner_product(vec1, vec2):
@@ -384,3 +392,53 @@ def vec(mat):
 
     """
     return matrix(mat, (len(mat), 1))
+
+
+def condition_number(mat):
+    """
+    Return the condition number of the given matrix.
+
+    The condition number of a matrix quantifies how hard it is to do
+    numerical computation with that matrix. It is usually defined as
+    the ratio of the norm of the matrix to the norm of its inverse, and
+    therefore depends on the norm used. One way to compute the condition
+    number with respect to the 2-norm is as the ratio of the matrix's
+    largest and smallest singular values. Since we have easy access to
+    those singular values, that is the algorithm we use.
+
+    The larger the condition number is, the worse the matrix is.
+
+    Parameters
+    ----------
+    mat : matrix
+        The matrix whose condition number you want.
+
+    Returns
+    -------
+
+    float
+        The nonnegative condition number of ``mat``.
+
+    Examples
+    --------
+
+    >>> condition_number(identity(3, typecode='d'))
+    1.0
+
+    >>> A = matrix([[2,1],[1,2]], tc='d')
+    >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
+    True
+
+    >>> A = matrix([[2,1j],[-1j,2]], tc='z')
+    >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
+    True
+
+    """
+    num_eigs = min(mat.size)
+    eigs = matrix(0, (num_eigs, 1), tc='d')
+    gesdd(mat, eigs)
+
+    if len(eigs) > 0:
+        return eigs[0]/eigs[-1]
+    else:
+        return 0