]> gitweb.michael.orlitzky.com - dunshire.git/blobdiff - dunshire/matrices.py
A bunch more doc fixes.
[dunshire.git] / dunshire / matrices.py
index bcf83778d62752436f7894b5a6fbad4cc3e1012e..f35827d16be600cec3ee5d1c2d0fbcea27a9484d 100644 (file)
@@ -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
@@ -17,8 +18,11 @@ def append_col(left, right):
     Parameters
     ----------
 
     Parameters
     ----------
 
-    left, right : matrix
-        The two matrices to append to one another.
+    left : matrix
+        The "original" matrix, the one that will wind up on the left.
+
+    right : matrix
+        The matrix to be appended on the right of ``left``.
 
     Returns
     -------
 
     Returns
     -------
@@ -56,8 +60,11 @@ def append_row(top, bottom):
     Parameters
     ----------
 
     Parameters
     ----------
 
-    top, bottom : matrix
-        The two matrices to append to one another.
+    top : matrix
+        The "original" matrix, the one that will wind up on top.
+
+    bottom : matrix
+        The matrix to be appended below ``top``.
 
     Returns
     -------
 
     Returns
     -------
@@ -142,7 +149,7 @@ def eigenvalues(symmat):
     eigs = matrix(0, (domain_dim, 1), tc='d')
 
     # Create a copy of ``symmat`` here because ``syevr`` clobbers it.
     eigs = matrix(0, (domain_dim, 1), tc='d')
 
     # Create a copy of ``symmat`` here because ``syevr`` clobbers it.
-    dummy = matrix(symmat, symmat.size)
+    dummy = copy(symmat)
     syevr(dummy, eigs)
     return list(eigs)
 
     syevr(dummy, eigs)
     return list(eigs)
 
@@ -264,8 +271,11 @@ def inner_product(vec1, vec2):
     Parameters
     ----------
 
     Parameters
     ----------
 
-    vec1, vec2 : matrix
-        The two vectors whose inner product you want.
+    vec1 : matrix
+        The first vector, whose inner product with ``vec2`` you want.
+
+    vec2 : matrix
+        The second vector, whose inner product with ``vec1`` you want.
 
     Returns
     -------
 
     Returns
     -------
@@ -329,17 +339,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 Euclidean 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``.
@@ -421,25 +467,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)
     eigs = matrix(0, (num_eigs, 1), tc='d')
     >>> 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)
+    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]