Add the matrix_of_transformation function to matrix_vector.py.

diff --git a/mjo/matrix_vector.py b/mjo/matrix_vector.py
index e1d1931e6c574053827038a9764bebf9ef0c482e..49d112a9f832e8849428abab5955a9f779d884d2 100644 (file)
--- a/mjo/matrix_vector.py
+++ b/mjo/matrix_vector.py
@@ -48,3 +48,107 @@ def isomorphism(matrix_space):
         return matrix_space(v.list())
 
     return (phi, phi_inverse)
+
+
+
+def matrix_of_transformation(T, V):
+    """
+    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.
+
+    INPUT:
+
+    - ``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``.
+
+    - ``V`` -- The vector or matrix space on which ``T`` is defined.
+
+    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``.
+
+    EXAMPLES:
+
+    The matrix of a transformation on a simple vector space should be
+    the expected matrix::
+
+        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]
+
+    A more complicated example confirms that we get a matrix consistent
+    with our ``matrix_to_vector`` function::
+
+        sage: M = MatrixSpace(QQ,3,3)
+        sage: Q = M([[0,1,0],[1,0,0],[0,0,1]])
+        sage: def f(x):
+        ....:     return Q*x*Q.inverse()
+        ....:
+        sage: F = matrix_of_transformation(f, M)
+        sage: F
+        [0 0 0 0 1 0 0 0 0]
+        [0 0 0 1 0 0 0 0 0]
+        [0 0 0 0 0 1 0 0 0]
+        [0 1 0 0 0 0 0 0 0]
+        [1 0 0 0 0 0 0 0 0]
+        [0 0 1 0 0 0 0 0 0]
+        [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: 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 = 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(0,n):
+        for i in range(0,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