Add the specnorm() function to the matrices module.

[dunshire.git] / dunshire / matrices.py

diff --git a/dunshire/matrices.py b/dunshire/matrices.py
index 2d4bb17c98a9187e2cacaa0094ef9c0cfe3c4600..a5a02e9d069eea87b2cf44cb70054133f3e4ae8c 100644 (file)
--- a/dunshire/matrices.py
+++ b/dunshire/matrices.py
@@ -3,6 +3,7 @@ Utility functions for working with CVXOPT matrices (instances of the
 class:`cvxopt.base.matrix` class).
 """
 
 class:`cvxopt.base.matrix` class).
 """
 

+from copy import copy

 from math import sqrt
 from cvxopt import matrix
 from cvxopt.lapack import gees, gesdd, syevr
 from math import sqrt
 from cvxopt import matrix
 from cvxopt.lapack import gees, gesdd, syevr


@@ -140,7 +141,10 @@ def eigenvalues(symmat):
 
     domain_dim = symmat.size[0]
     eigs = matrix(0, (domain_dim, 1), tc='d')
 
     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)
 
 
     return list(eigs)
 
 


@@ -326,17 +330,53 @@ def norm(matrix_or_vector):
     --------
 
         >>> v = matrix([1,1])
     --------
 
         >>> v = matrix([1,1])

-        >>> print('{:.5f}'.format(norm(v)))
-        1.41421
+        >>> norm(v)
+        1.414...

 
         >>> A = matrix([1,1,1,1], (2,2))
         >>> norm(A)
 
         >>> A = matrix([1,1,1,1], (2,2))
         >>> norm(A)

-        2.0
+        2.0...

 
     """
     return sqrt(inner_product(matrix_or_vector, matrix_or_vector))
 
 
 
     """
     return sqrt(inner_product(matrix_or_vector, matrix_or_vector))
 
 

+def specnorm(mat):
+    """
+    Return the spectral norm of a matrix.
+
+    The spectral norm of a matrix is its largest singular value, and it
+    corresponds to the operator norm induced by the vector ``2``-norm.
+
+    Parameters
+    ----------
+
+    mat : matrix
+        The matrix whose spectral norm you want.
+
+    Examples:
+
+        >>> specnorm(identity(3))
+        1.0
+
+        >>> specnorm(5*identity(4))
+        5.0
+
+    """
+    num_eigs = min(mat.size)
+    eigs = matrix(0, (num_eigs, 1), tc='d')
+    typecode = 'd'
+    if any([isinstance(entry, complex) for entry in mat]):
+        typecode = 'z'
+    dummy = matrix(mat, mat.size, tc=typecode)
+    gesdd(dummy, eigs)
+
+    if len(eigs) > 0:
+        return eigs[0]
+    else:
+        return 0
+
+

 def vec(mat):
     """
     Create a long vector in column-major order from ``mat``.
 def vec(mat):
     """
     Create a long vector in column-major order from ``mat``.


@@ -418,25 +458,25 @@ def condition_number(mat):
     Examples
     --------
 
     Examples
     --------
 

-    >>> condition_number(identity(1, typecode='d'))
-    1.0
-    >>> condition_number(identity(2, typecode='d'))
-    1.0
-    >>> condition_number(identity(3, typecode='d'))
+    >>> condition_number(identity(3))

     1.0
 
     1.0
 

-    >>> A = matrix([[2,1],[1,2]], tc='d')
+    >>> A = matrix([[2,1],[1,2]])

     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
     True
 
     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
     True
 

-    >>> A = matrix([[2,1j],[-1j,2]], tc='z')
+    >>> A = matrix([[2,1j],[-1j,2]])

     >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
     True
 
     """
     num_eigs = min(mat.size)
     >>> 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)
+    eigs = matrix(0, (num_eigs, 1), tc='d')
+    typecode = 'd'
+    if any([isinstance(entry, complex) for entry in mat]):
+        typecode = 'z'
+    dummy = matrix(mat, mat.size, tc=typecode)
+    gesdd(dummy, eigs)

 
     if len(eigs) > 0:
         return eigs[0]/eigs[-1]
 
     if len(eigs) > 0:
         return eigs[0]/eigs[-1]