mjo/matrix_vector: replace matrix_of_transformation with basis_repr_of_operator.

diff --git a/mjo/matrix_vector.py b/mjo/matrix_vector.py
index 62b75fd43e2299d605d7000f7c741af6e1bf6690..83bf313a82cde7c689d10337efff8054301acad7 100644 (file)
--- a/mjo/matrix_vector.py
+++ b/mjo/matrix_vector.py
@@ -127,45 +127,59 @@ def basis_representation(M):
 
 
 
-def matrix_of_transformation(T, V):
+def basis_repr_of_operator(M, L):
     """
-    Compute the matrix of a linear transformation ``T``, `$T : V
-    \rightarrow V$` with domain/range ``V``. This essentially uses the
-    Riesz representation theorem to represent the entries of the matrix
-    of ``T`` in terms of inner products.
+    Return the matrix of the operator `L` with respect to the basis
+    `M` if `M` is a list of basis vectors for a matrix space; or with
+    respect to the standard basis of `M` if `M` is a ``MatrixSpace``.
 
-    INPUT:
+    This function is necessary because SageMath does not know that
+    matrix spaces are vector spaces, and it moreover it doesn't know
+    that (for example) the subspace of symmetric matrices is a matrix
+    space in its own right.
 
-    - ``T`` -- The linear transformation whose matrix we should
-               compute. This should be a callable function that
-               takes as its single argument an element of ``V``.
+    Use ``linear_transformation().matrix()`` instead if you have a
+    true ``VectorSpace``.
 
-    - ``V`` -- The vector or matrix space on which ``T`` is defined.
+    INPUT:
+
+    - ``M`` -- Either a ``MatrixSpace``, or a list of matrices that form
+               a basis for a matrix space.
 
     OUTPUT:
 
-    If the dimension of ``V`` is `$n$`, we return an `$n \times n$`
-    matrix that represents ``T`` with respect to the standard basis of
-    ``V``.
+    If the matrix space associated with `M` has dimension `n`, then an
+    `n`-by-`n` matrix over the same base ring is returned.
 
     SETUP::
 
         sage: from mjo.matrix_vector import (basis_representation,
-        ....:                                matrix_of_transformation)
+        ....:                                basis_repr_of_operator)
 
     EXAMPLES:
 
-    The matrix of a transformation on a simple vector space should be
-    the expected matrix::
+    The matrix of the identity operator on the space of two-by-two
+    symmetric matrices is the identity matrix, regardless of the basis::
 
-        sage: V = VectorSpace(QQ, 3)
-        sage: def f(x):
-        ....:     return 3*x
-        ....:
-        sage: matrix_of_transformation(f, V)
-        [3 0 0]
-        [0 3 0]
-        [0 0 3]
+        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: identity = lambda X: X
+        sage: basis_repr_of_operator(basis, identity)
+        [1 0 0]
+        [0 1 0]
+        [0 0 1]
+        sage: E11 = matrix(QQ,[[2,0],[0,0]])
+        sage: E12 = matrix(QQ,[[0,2],[2,0]])
+        sage: basis = [E11, E12, E22]
+        sage: basis_repr_of_operator(basis, identity)
+        [1 0 0]
+        [0 1 0]
+        [0 0 1]
 
     A more complicated example confirms that we get a matrix consistent
     with our ``matrix_to_vector`` function::
@@ -175,7 +189,7 @@ def matrix_of_transformation(T, V):
         sage: def f(x):
         ....:     return Q*x*Q.inverse()
         ....:
-        sage: F = matrix_of_transformation(f, M)
+        sage: F = basis_repr_of_operator(M, f)
         sage: F
         [0 0 0 0 1 0 0 0 0]
         [0 0 0 1 0 0 0 0 0]
@@ -186,50 +200,27 @@ def matrix_of_transformation(T, V):
         [0 0 0 0 0 0 0 1 0]
         [0 0 0 0 0 0 1 0 0]
         [0 0 0 0 0 0 0 0 1]
-        sage: phi, phi_inv = isomorphism(M)
+        sage: phi, phi_inv = basis_representation(M)
         sage: X = M([[1,2,3],[4,5,6],[7,8,9]])
-        sage: F*phi(X)
-        (5, 4, 6, 2, 1, 3, 8, 7, 9)
-        sage: phi(f(X))
-        (5, 4, 6, 2, 1, 3, 8, 7, 9)
         sage: F*phi(X) == phi(f(X))
         True
 
     """
-    n = V.dimension()
-    B = list(V.basis())
-
-    def inner_product(v, w):
-        # An inner product function that works for both matrices and
-        # vectors.
-        if callable(getattr(v, 'inner_product', None)):
-            return v.inner_product(w)
-        elif callable(getattr(v, 'matrix_space', None)):
-            # V must be a matrix space?
-            return (v*w.transpose()).trace()
-        else:
-            raise ValueError('inner_product only works on vectors and matrices')
-
-    def apply(L, x):
-        # A "call" that works for both matrices and functions.
-        if callable(getattr(L, 'matrix_space', None)):
-            # L is a matrix, and we need to use "multiply" to call it.
-            return L*x
-        else:
-            # If L isn't a matrix, try this. It works for python
-            # functions at least.
-            return L(x)
-
-    entries = []
-    for j in range(n):
-        for i in range(n):
-            entry = inner_product(apply(T,B[i]), B[j])
-            entries.append(entry)
-
-    # Construct the matrix space in which our return value will lie.
-    W = MatrixSpace(V.base_ring(), n, n)
-
-    # And make a matrix out of our list of entries.
-    A = W(entries)
-
-    return A
+    if is_MatrixSpace(M):
+        basis_space = M
+        basis = list(M.basis())
+    else:
+        basis_space = M[0].matrix_space()
+        basis = M
+
+    (phi, phi_inv) = basis_representation(M)
+
+    # Get a basis for the image space. Since phi is an isometry,
+    # it takes one basis to another.
+    image_basis = [ phi(b) for b in basis ]
+
+    # Now construct the image space itself equipped with our custom basis.
+    W = VectorSpace(basis_space.base_ring(), len(basis))
+    W = W.span_of_basis(image_basis)
+
+    return matrix.column( W.coordinates(phi(L(b))) for b in basis )