mjo/matrix_vector.py

   1 """
   2 There is an explicit isomorphism between all finite-dimensional vector
   3 spaces. In particular, there is an isomorphism between the m-by-n
   4 matrices and `$R^(m \times n)$`. Since most vector operations are not
   5 available on Sage matrices, we have to go back and forth between these
   6 two vector spaces often.
   7 """
   8
   9 from sage.all import *
  10 from sage.matrix.matrix_space import is_MatrixSpace
  11
  12 def _mat2vec(m):
  13     return vector(m.base_ring(), m.list())
  14
  15 def basis_representation(M):
  16     """
  17     Return the forward (``MatrixSpace`` -> ``VectorSpace``) and
  18     inverse isometries, as a pair, that take elements of the given
  19     ``MatrixSpace`` `M` to their representations as "long vectors,"
  20     and vice-versa.
  21
  22     The argument ``M`` can be either a ``MatrixSpace`` or a basis for
  23     a space of matrices. This function is needed because SageMath does
  24     not know that matrix spaces are vector spaces, and therefore
  25     cannot perform common operations with them -- like computing the
  26     basis representation of an element.
  27
  28     Moreover, the ability to pass in a basis (rather than a
  29     ``MatrixSpace``) is needed because SageMath has no way to express
  30     that e.g. a (sub)space of symmetric matrices is itself a
  31     ``MatrixSpace``.
  32
  33     INPUT:
  34
  35     - ``M`` -- Either a ``MatrixSpace``, or a list of matrices that form
  36                a basis for a matrix space.
  37
  38     OUTPUT:
  39
  40     A pair of isometries ``(phi, phi_inv)``.
  41
  42     If the matrix space associated with `M` has dimension `n`, then
  43     ``phi`` will map its elements to vectors of length `n` over the
  44     same base ring. The inverse map ``phi_inv`` reverses that
  45     operation.
  46
  47     SETUP::
  48
  49         sage: from mjo.matrix_vector import basis_representation
  50
  51     EXAMPLES:
  52
  53     This function computes the correct coordinate representations (of
  54     length 3) for a basis of the space of two-by-two symmetric
  55     matrices, the the inverse does indeed invert the process::
  56
  57         sage: E11 = matrix(QQbar,[ [1,0],
  58         ....:                      [0,0] ])
  59         sage: E12 = matrix(QQbar,[ [0, 1/sqrt(2)],
  60         ....:                      [1/sqrt(2), 0] ])
  61         sage: E22 = matrix(QQbar,[ [0,0],
  62         ....:                      [0,1] ])
  63         sage: basis = [E11, E12, E22]
  64         sage: phi, phi_inv = basis_representation(basis)
  65         sage: phi(E11); phi(E12); phi(E22)
  66         (1, 0, 0)
  67         (0, 1, 0)
  68         (0, 0, 1)
  69         sage: phi_inv(phi(E11)) == E11
  70         True
  71         sage: phi_inv(phi(E12)) == E12
  72         True
  73         sage: phi_inv(phi(E22)) == E22
  74         True
  75
  76     MatrixSpace arguments work too::
  77
  78         sage: M = MatrixSpace(QQ,2)
  79         sage: phi, phi_inv = basis_representation(M)
  80         sage: X = matrix(QQ, [ [1,2],
  81         ....:                  [3,4] ])
  82         sage: phi(X)
  83         (1, 2, 3, 4)
  84         sage: phi_inv(phi(X)) == X
  85         True
  86
  87     TESTS:
  88
  89     The inverse is generally an inverse::
  90
  91         sage: set_random_seed()
  92         sage: n = ZZ.random_element(10)
  93         sage: M = MatrixSpace(QQ,n)
  94         sage: X = M.random_element()
  95         sage: (phi, phi_inv) = basis_representation(M)
  96         sage: phi_inv(phi(X)) == X
  97         True
  98
  99     """
 100     if is_MatrixSpace(M):
 101         basis_space = M
 102         basis = list(M.basis())
 103     else:
 104         basis_space = M[0].matrix_space()
 105         basis = M
 106
 107     def phi(X):
 108         """
 109         The isometry sending ``X`` to its representation as a long vector.
 110         """
 111         if X not in basis_space:
 112             raise ValueError("X does not live in the domain of phi")
 113
 114         V = VectorSpace(basis_space.base_ring(), X.nrows()*X.ncols())
 115         W = V.span_of_basis( _mat2vec(s) for s in basis )
 116         return W.coordinate_vector(_mat2vec(X))
 117
 118     def phi_inv(Y):
 119         """
 120         The isometry sending the long vector `Y` to an element of either
 121         `M` or the span of `M` (depending on whether or not ``M``
 122         is a ``MatrixSpace`` or a basis).
 123         """
 124         return sum( Y[i]*basis[i] for i in range(len(Y)) )
 125
 126     return (phi, phi_inv)
 127
 128
 129
 130 def matrix_of_transformation(T, V):
 131     """
 132     Compute the matrix of a linear transformation ``T``, `$T : V
 133     \rightarrow V$` with domain/range ``V``. This essentially uses the
 134     Riesz representation theorem to represent the entries of the matrix
 135     of ``T`` in terms of inner products.
 136
 137     INPUT:
 138
 139     - ``T`` -- The linear transformation whose matrix we should
 140                compute. This should be a callable function that
 141                takes as its single argument an element of ``V``.
 142
 143     - ``V`` -- The vector or matrix space on which ``T`` is defined.
 144
 145     OUTPUT:
 146
 147     If the dimension of ``V`` is `$n$`, we return an `$n \times n$`
 148     matrix that represents ``T`` with respect to the standard basis of
 149     ``V``.
 150
 151     SETUP::
 152
 153         sage: from mjo.matrix_vector import (basis_representation,
 154         ....:                                matrix_of_transformation)
 155
 156     EXAMPLES:
 157
 158     The matrix of a transformation on a simple vector space should be
 159     the expected matrix::
 160
 161         sage: V = VectorSpace(QQ, 3)
 162         sage: def f(x):
 163         ....:     return 3*x
 164         ....:
 165         sage: matrix_of_transformation(f, V)
 166         [3 0 0]
 167         [0 3 0]
 168         [0 0 3]
 169
 170     A more complicated example confirms that we get a matrix consistent
 171     with our ``matrix_to_vector`` function::
 172
 173         sage: M = MatrixSpace(QQ,3,3)
 174         sage: Q = M([[0,1,0],[1,0,0],[0,0,1]])
 175         sage: def f(x):
 176         ....:     return Q*x*Q.inverse()
 177         ....:
 178         sage: F = matrix_of_transformation(f, M)
 179         sage: F
 180         [0 0 0 0 1 0 0 0 0]
 181         [0 0 0 1 0 0 0 0 0]
 182         [0 0 0 0 0 1 0 0 0]
 183         [0 1 0 0 0 0 0 0 0]
 184         [1 0 0 0 0 0 0 0 0]
 185         [0 0 1 0 0 0 0 0 0]
 186         [0 0 0 0 0 0 0 1 0]
 187         [0 0 0 0 0 0 1 0 0]
 188         [0 0 0 0 0 0 0 0 1]
 189         sage: phi, phi_inv = isomorphism(M)
 190         sage: X = M([[1,2,3],[4,5,6],[7,8,9]])
 191         sage: F*phi(X)
 192         (5, 4, 6, 2, 1, 3, 8, 7, 9)
 193         sage: phi(f(X))
 194         (5, 4, 6, 2, 1, 3, 8, 7, 9)
 195         sage: F*phi(X) == phi(f(X))
 196         True
 197
 198     """
 199     n = V.dimension()
 200     B = list(V.basis())
 201
 202     def inner_product(v, w):
 203         # An inner product function that works for both matrices and
 204         # vectors.
 205         if callable(getattr(v, 'inner_product', None)):
 206             return v.inner_product(w)
 207         elif callable(getattr(v, 'matrix_space', None)):
 208             # V must be a matrix space?
 209             return (v*w.transpose()).trace()
 210         else:
 211             raise ValueError('inner_product only works on vectors and matrices')
 212
 213     def apply(L, x):
 214         # A "call" that works for both matrices and functions.
 215         if callable(getattr(L, 'matrix_space', None)):
 216             # L is a matrix, and we need to use "multiply" to call it.
 217             return L*x
 218         else:
 219             # If L isn't a matrix, try this. It works for python
 220             # functions at least.
 221             return L(x)
 222
 223     entries = []
 224     for j in range(n):
 225         for i in range(n):
 226             entry = inner_product(apply(T,B[i]), B[j])
 227             entries.append(entry)
 228
 229     # Construct the matrix space in which our return value will lie.
 230     W = MatrixSpace(V.base_ring(), n, n)
 231
 232     # And make a matrix out of our list of entries.
 233     A = W(entries)
 234
 235     return A