mjo/basis_repr.py

   1 r"""
   2 In an `n`-dimensional vector space, representation with respect to
   3 a basis is an isometry between that space and `\mathbb{R}^{n}`.
   4
   5 Sage is able to go back/forth relatively easy when you start with a
   6 ``VectorSpace``, but unfortunately, it does not know that a
   7 ``MatrixSpace`` is also a ``VectorSpace``. So, this module exists to
   8 perform the "basis representation" isometry between a matrix space and
   9 a vector space of the same dimension.
  10
  11 """
  12
  13 from sage.all import *
  14 from sage.matrix.matrix_space import is_MatrixSpace
  15
  16 def _mat2vec(m):
  17     return vector(m.base_ring(), m.list())
  18
  19 def basis_repr(M):
  20     """
  21     Return the forward (``MatrixSpace`` -> ``VectorSpace``) and
  22     inverse isometries, as a pair, that take elements of the given
  23     ``MatrixSpace`` `M` to their representations as "long vectors,"
  24     and vice-versa.
  25
  26     The argument ``M`` can be either a ``MatrixSpace`` or a basis for
  27     a space of matrices. This function is needed because SageMath does
  28     not know that matrix spaces are vector spaces, and therefore
  29     cannot perform common operations with them -- like computing the
  30     basis representation of an element.
  31
  32     Moreover, the ability to pass in a basis (rather than a
  33     ``MatrixSpace``) is needed because SageMath has no way to express
  34     that e.g. a (sub)space of symmetric matrices is itself a
  35     ``MatrixSpace``.
  36
  37     INPUT:
  38
  39     - ``M`` -- Either a ``MatrixSpace``, or a list of matrices that form
  40                a basis for a matrix space.
  41
  42     OUTPUT:
  43
  44     A pair of isometries ``(phi, phi_inv)``.
  45
  46     If the matrix space associated with `M` has dimension `n`, then
  47     ``phi`` will map its elements to vectors of length `n` over the
  48     same base ring. The inverse map ``phi_inv`` reverses that
  49     operation.
  50
  51     SETUP::
  52
  53         sage: from mjo.basis_repr import basis_repr
  54
  55     EXAMPLES:
  56
  57     This function computes the correct coordinate representations (of
  58     length 3) for a basis of the space of two-by-two symmetric
  59     matrices, the the inverse does indeed invert the process::
  60
  61         sage: E11 = matrix(QQbar,[ [1,0],
  62         ....:                      [0,0] ])
  63         sage: E12 = matrix(QQbar,[ [0, 1/sqrt(2)],
  64         ....:                      [1/sqrt(2), 0] ])
  65         sage: E22 = matrix(QQbar,[ [0,0],
  66         ....:                      [0,1] ])
  67         sage: basis = [E11, E12, E22]
  68         sage: phi, phi_inv = basis_repr(basis)
  69         sage: phi(E11); phi(E12); phi(E22)
  70         (1, 0, 0)
  71         (0, 1, 0)
  72         (0, 0, 1)
  73         sage: phi_inv(phi(E11)) == E11
  74         True
  75         sage: phi_inv(phi(E12)) == E12
  76         True
  77         sage: phi_inv(phi(E22)) == E22
  78         True
  79
  80     MatrixSpace arguments work too::
  81
  82         sage: M = MatrixSpace(QQ,2)
  83         sage: phi, phi_inv = basis_repr(M)
  84         sage: X = matrix(QQ, [ [1,2],
  85         ....:                  [3,4] ])
  86         sage: phi(X)
  87         (1, 2, 3, 4)
  88         sage: phi_inv(phi(X)) == X
  89         True
  90
  91     TESTS:
  92
  93     The inverse is generally an inverse::
  94
  95         sage: n = ZZ.random_element(10)
  96         sage: M = MatrixSpace(QQ,n)
  97         sage: X = M.random_element()
  98         sage: (phi, phi_inv) = basis_repr(M)
  99         sage: phi_inv(phi(X)) == X
 100         True
 101
 102     """
 103     if is_MatrixSpace(M):
 104         basis_space = M
 105         basis = list(M.basis())
 106     else:
 107         basis_space = M[0].matrix_space()
 108         basis = M
 109
 110     def phi(X):
 111         """
 112         The isometry sending ``X`` to its representation as a long vector.
 113         """
 114         if X not in basis_space:
 115             raise ValueError("X does not live in the domain of phi")
 116
 117         V = VectorSpace(basis_space.base_ring(), X.nrows()*X.ncols())
 118         W = V.span_of_basis( _mat2vec(s) for s in basis )
 119         return W.coordinate_vector(_mat2vec(X))
 120
 121     def phi_inv(Y):
 122         """
 123         The isometry sending the long vector `Y` to an element of either
 124         `M` or the span of `M` (depending on whether or not ``M``
 125         is a ``MatrixSpace`` or a basis).
 126         """
 127         return basis_space.linear_combination( zip(Y,basis) )
 128
 129     return (phi, phi_inv)
 130
 131
 132
 133 def basis_repr_of_operator(M, L):
 134     """
 135     Return the matrix of the operator `L` with respect to the basis
 136     `M` if `M` is a list of basis vectors for a matrix space; or with
 137     respect to the standard basis of `M` if `M` is a ``MatrixSpace``.
 138
 139     This function is necessary because SageMath does not know that
 140     matrix spaces are vector spaces, and it moreover it doesn't know
 141     that (for example) the subspace of symmetric matrices is a matrix
 142     space in its own right.
 143
 144     Use ``linear_transformation().matrix()`` instead if you have a
 145     true ``VectorSpace``.
 146
 147     INPUT:
 148
 149     - ``M`` -- Either a ``MatrixSpace``, or a list of matrices that form
 150                a basis for a matrix space.
 151
 152     OUTPUT:
 153
 154     If the matrix space associated with `M` has dimension `n`, then an
 155     `n`-by-`n` matrix over the same base ring is returned.
 156
 157     SETUP::
 158
 159         sage: from mjo.basis_repr import (basis_repr,
 160         ....:                             basis_repr_of_operator)
 161
 162     EXAMPLES:
 163
 164     The matrix of the identity operator on the space of two-by-two
 165     symmetric matrices is the identity matrix, regardless of the basis::
 166
 167         sage: E11 = matrix(QQbar,[ [1,0],
 168         ....:                      [0,0] ])
 169         sage: E12 = matrix(QQbar,[ [0, 1/sqrt(2)],
 170         ....:                      [1/sqrt(2), 0] ])
 171         sage: E22 = matrix(QQbar,[ [0,0],
 172         ....:                      [0,1] ])
 173         sage: basis = [E11, E12, E22]
 174         sage: identity = lambda X: X
 175         sage: basis_repr_of_operator(basis, identity)
 176         [1 0 0]
 177         [0 1 0]
 178         [0 0 1]
 179         sage: E11 = matrix(QQ,[[2,0],[0,0]])
 180         sage: E12 = matrix(QQ,[[0,2],[2,0]])
 181         sage: basis = [E11, E12, E22]
 182         sage: basis_repr_of_operator(basis, identity)
 183         [1 0 0]
 184         [0 1 0]
 185         [0 0 1]
 186
 187     A more complicated example confirms that we get a matrix consistent
 188     with our ``matrix_to_vector`` function::
 189
 190         sage: M = MatrixSpace(QQ,3,3)
 191         sage: Q = M([[0,1,0],[1,0,0],[0,0,1]])
 192         sage: def f(x):
 193         ....:     return Q*x*Q.inverse()
 194         ....:
 195         sage: F = basis_repr_of_operator(M, f)
 196         sage: F
 197         [0 0 0 0 1 0 0 0 0]
 198         [0 0 0 1 0 0 0 0 0]
 199         [0 0 0 0 0 1 0 0 0]
 200         [0 1 0 0 0 0 0 0 0]
 201         [1 0 0 0 0 0 0 0 0]
 202         [0 0 1 0 0 0 0 0 0]
 203         [0 0 0 0 0 0 0 1 0]
 204         [0 0 0 0 0 0 1 0 0]
 205         [0 0 0 0 0 0 0 0 1]
 206         sage: phi, phi_inv = basis_repr(M)
 207         sage: X = M([[1,2,3],[4,5,6],[7,8,9]])
 208         sage: F*phi(X) == phi(f(X))
 209         True
 210
 211     """
 212     if is_MatrixSpace(M):
 213         basis_space = M
 214         basis = list(M.basis())
 215     else:
 216         basis_space = M[0].matrix_space()
 217         basis = M
 218
 219     (phi, phi_inv) = basis_repr(M)
 220
 221     # Get a basis for the image space. Since phi is an isometry,
 222     # it takes one basis to another.
 223     image_basis = [ phi(b) for b in basis ]
 224
 225     # Now construct the image space itself equipped with our custom basis.
 226     W = VectorSpace(basis_space.base_ring(), len(basis))
 227     W = W.span_of_basis(image_basis)
 228
 229     return matrix.column( W.coordinates(phi(L(b))) for b in basis )