mjo/matrix_vector: replace isomorphism with basis_representation().

diff --git a/mjo/cone/doubly_nonnegative.py b/mjo/cone/doubly_nonnegative.py
index 1fb350cb46559a55b982ee3de15678e02d5a04de..373557c6c72dde4974dfdce30efe435d0d526cc2 100644 (file)
--- a/mjo/cone/doubly_nonnegative.py
+++ b/mjo/cone/doubly_nonnegative.py
@@ -16,7 +16,7 @@ from sage.all import *
 from mjo.cone.symmetric_psd import (factor_psd,
                                     is_symmetric_psd,
                                     random_symmetric_psd)
-from mjo.matrix_vector import isomorphism
+from mjo.matrix_vector import basis_representation
 
 
 def is_doubly_nonnegative(A):
@@ -383,7 +383,7 @@ def is_extreme_doubly_nonnegative(A):
     # can't compute the dimension of a set of matrices anyway, so we
     # convert them all to vectors and just ask for the dimension of the
     # resulting vector space.
-    (phi, phi_inverse) = isomorphism(A.matrix_space())
+    (phi, phi_inverse) = basis_representation(A.matrix_space())
     vectors = map(phi,spanning_set)
 
     V = span(vectors, A.base_ring())

diff --git a/mjo/matrix_vector.py b/mjo/matrix_vector.py
index 111712b863d3a0dcf9f730798a5efc9fcd236de6..4b4818df133a54fc16c6cf7e22b9a0e8eda87796 100644 (file)
--- a/mjo/matrix_vector.py
+++ b/mjo/matrix_vector.py
@@ -7,51 +7,123 @@ two vector spaces often.
 """
 
 from sage.all import *
+from sage.matrix.matrix_space import is_MatrixSpace
 
-def isomorphism(matrix_space):
+def _mat2vec(m):
+    return vector(m.base_ring(), m.list())
+
+def basis_representation(M):
     """
-    Create isomorphism (i.e. the function) that converts elements
-    of a matrix space into those of the corresponding finite-dimensional
-    vector space.
+    Return the forward (``MatrixSpace`` -> ``VectorSpace``) and
+    inverse isometries, as a pair, that take elements of the given
+    ``MatrixSpace`` `M` to their representations as "long vectors,"
+    and vice-versa.
+
+    The argument ``M`` can be either a ``MatrixSpace`` or a basis for
+    a space of matrices. This function is needed because SageMath does
+    not know that matrix spaces are vector spaces, and therefore
+    cannot perform common operations with them -- like computing the
+    basis representation of an element.
+
+    Moreover, the ability to pass in a basis (rather than a
+    ``MatrixSpace``) is needed because SageMath has no way to express
+    that e.g. a (sub)space of symmetric matrices is itself a
+    ``MatrixSpace``.
 
     INPUT:
 
-    - matrix_space: A finite-dimensional ``MatrixSpace`` object.
+    - ``M`` -- Either a ``MatrixSpace``, or a list of matrices that form
+               a basis for a matrix space.
 
     OUTPUT:
 
-    - (phi, phi_inverse): If ``matrix_space`` has dimension m*n, then
-                          ``phi`` will map m-by-n matrices to R^(m*n).
-                          The inverse mapping ``phi_inverse`` will go
-                          the other way.
+    A pair of isometries ``(phi, phi_inv)``.
+
+    If the matrix space associated with `M` has dimension `n`, then
+    ``phi`` will map its elements to vectors of length `n` over the
+    same base ring. The inverse map ``phi_inv`` reverses that
+    operation.
 
     SETUP::
 
-        sage: from mjo.matrix_vector import isomorphism
+        sage: from mjo.matrix_vector import basis_representation
+
+    EXAMPLES:
+
+    This function computes the correct coordinate representations (of
+    length 3) for a basis of the space of two-by-two symmetric
+    matrices, the the inverse does indeed invert the process::
+
+        sage: E11 = matrix(QQbar,[ [1,0],
+        ....:                      [0,0] ])
+        sage: E12 = matrix(QQbar,[ [0, 1/sqrt(2)],
+        ....:                      [1/sqrt(2), 0] ])
+        sage: E22 = matrix(QQbar,[ [0,0],
+        ....:                      [0,1] ])
+        sage: basis = [E11, E12, E22]
+        sage: phi, phi_inv = basis_representation(basis)
+        sage: phi(E11); phi(E12); phi(E22)
+        (1, 0, 0)
+        (0, 1, 0)
+        (0, 0, 1)
+        sage: phi_inv(phi(E11)) == E11
+        True
+        sage: phi_inv(phi(E12)) == E12
+        True
+        sage: phi_inv(phi(E22)) == E22
+        True
+
+    MatrixSpace arguments work too::
+
+        sage: M = MatrixSpace(QQ,2)
+        sage: phi, phi_inv = basis_representation(M)
+        sage: X = matrix(QQ, [ [1,2],
+        ....:                  [3,4] ])
+        sage: phi(X)
+        (1, 2, 3, 4)
+        sage: phi_inv(phi(X)) == X
+        True
+
+    TESTS:
 
-    EXAMPLES::
+    The inverse is generally an inverse::
 
-        sage: M = MatrixSpace(QQ,4,4)
-        sage: (p, p_inv) = isomorphism(M)
-        sage: m = M(range(16))
-        sage: p_inv(p(m)) == m
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(10)
+        sage: M = MatrixSpace(QQ,n)
+        sage: X = M.random_element()
+        sage: (phi, phi_inv) = basis_representation(M)
+        sage: phi_inv(phi(X)) == X
         True
 
     """
-    from sage.matrix.matrix_space import is_MatrixSpace
-    if not is_MatrixSpace(matrix_space):
-        raise TypeError('argument must be a matrix space')
+    if is_MatrixSpace(M):
+        basis_space = M
+        basis = list(M.basis())
+    else:
+        basis_space = M[0].matrix_space()
+        basis = M
 
-    base_ring = matrix_space.base_ring()
-    vector_space = VectorSpace(base_ring, matrix_space.dimension())
+    def phi(X):
+        """
+        The isometry sending ``X`` to its representation as a long vector.
+        """
+        if X not in basis_space:
+            raise ValueError("X does not live in the domain of phi")
 
-    def phi(m):
-        return vector_space(m.list())
+        V = VectorSpace(basis_space.base_ring(), X.nrows()*X.ncols())
+        W = V.span_of_basis( _mat2vec(s) for s in basis )
+        return W.coordinate_vector(_mat2vec(X))
 
-    def phi_inverse(v):
-        return matrix_space(v.list())
+    def phi_inv(Y):
+        """
+        The isometry sending the long vector `Y` to an element of either
+        `M` or the span of `M` (depending on whether or not ``M``
+        is a ``MatrixSpace`` or a basis).
+        """
+        return sum( Y[i]*basis[i] for i in range(len(Y)) )
 
-    return (phi, phi_inverse)
+    return (phi, phi_inv)
 
 
 
@@ -78,7 +150,8 @@ def matrix_of_transformation(T, V):
 
     SETUP::
 
-        sage: from mjo.matrix_vector import isomorphism, matrix_of_transformation
+        sage: from mjo.matrix_vector import (basis_representation,
+        ....:                                matrix_of_transformation)
 
     EXAMPLES: