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: set_random_seed()
  96         sage: n = ZZ.random_element(10)
  97         sage: M = MatrixSpace(QQ,n)
  98         sage: X = M.random_element()
  99         sage: (phi, phi_inv) = basis_repr(M)
 100         sage: phi_inv(phi(X)) == X
 101         True
 102
 103     """
 104     if is_MatrixSpace(M):
 105         basis_space = M
 106         basis = list(M.basis())
 107     else:
 108         basis_space = M[0].matrix_space()
 109         basis = M
 110
 111     def phi(X):
 112         """
 113         The isometry sending ``X`` to its representation as a long vector.
 114         """
 115         if X not in basis_space:
 116             raise ValueError("X does not live in the domain of phi")
 117
 118         V = VectorSpace(basis_space.base_ring(), X.nrows()*X.ncols())
 119         W = V.span_of_basis( _mat2vec(s) for s in basis )
 120         return W.coordinate_vector(_mat2vec(X))
 121
 122     def phi_inv(Y):
 123         """
 124         The isometry sending the long vector `Y` to an element of either
 125         `M` or the span of `M` (depending on whether or not ``M``
 126         is a ``MatrixSpace`` or a basis).
 127         """
 128         return basis_space.linear_combination( zip(Y,basis) )
 129
 130     return (phi, phi_inv)
 131
 132
 133
 134 def basis_repr_of_operator(M, L):
 135     """
 136     Return the matrix of the operator `L` with respect to the basis
 137     `M` if `M` is a list of basis vectors for a matrix space; or with
 138     respect to the standard basis of `M` if `M` is a ``MatrixSpace``.
 139
 140     This function is necessary because SageMath does not know that
 141     matrix spaces are vector spaces, and it moreover it doesn't know
 142     that (for example) the subspace of symmetric matrices is a matrix
 143     space in its own right.
 144
 145     Use ``linear_transformation().matrix()`` instead if you have a
 146     true ``VectorSpace``.
 147
 148     INPUT:
 149
 150     - ``M`` -- Either a ``MatrixSpace``, or a list of matrices that form
 151                a basis for a matrix space.
 152
 153     OUTPUT:
 154
 155     If the matrix space associated with `M` has dimension `n`, then an
 156     `n`-by-`n` matrix over the same base ring is returned.
 157
 158     SETUP::
 159
 160         sage: from mjo.basis_repr import (basis_repr,
 161         ....:                             basis_repr_of_operator)
 162
 163     EXAMPLES:
 164
 165     The matrix of the identity operator on the space of two-by-two
 166     symmetric matrices is the identity matrix, regardless of the basis::
 167
 168         sage: E11 = matrix(QQbar,[ [1,0],
 169         ....:                      [0,0] ])
 170         sage: E12 = matrix(QQbar,[ [0, 1/sqrt(2)],
 171         ....:                      [1/sqrt(2), 0] ])
 172         sage: E22 = matrix(QQbar,[ [0,0],
 173         ....:                      [0,1] ])
 174         sage: basis = [E11, E12, E22]
 175         sage: identity = lambda X: X
 176         sage: basis_repr_of_operator(basis, identity)
 177         [1 0 0]
 178         [0 1 0]
 179         [0 0 1]
 180         sage: E11 = matrix(QQ,[[2,0],[0,0]])
 181         sage: E12 = matrix(QQ,[[0,2],[2,0]])
 182         sage: basis = [E11, E12, E22]
 183         sage: basis_repr_of_operator(basis, identity)
 184         [1 0 0]
 185         [0 1 0]
 186         [0 0 1]
 187
 188     A more complicated example confirms that we get a matrix consistent
 189     with our ``matrix_to_vector`` function::
 190
 191         sage: M = MatrixSpace(QQ,3,3)
 192         sage: Q = M([[0,1,0],[1,0,0],[0,0,1]])
 193         sage: def f(x):
 194         ....:     return Q*x*Q.inverse()
 195         ....:
 196         sage: F = basis_repr_of_operator(M, f)
 197         sage: F
 198         [0 0 0 0 1 0 0 0 0]
 199         [0 0 0 1 0 0 0 0 0]
 200         [0 0 0 0 0 1 0 0 0]
 201         [0 1 0 0 0 0 0 0 0]
 202         [1 0 0 0 0 0 0 0 0]
 203         [0 0 1 0 0 0 0 0 0]
 204         [0 0 0 0 0 0 0 1 0]
 205         [0 0 0 0 0 0 1 0 0]
 206         [0 0 0 0 0 0 0 0 1]
 207         sage: phi, phi_inv = basis_repr(M)
 208         sage: X = M([[1,2,3],[4,5,6],[7,8,9]])
 209         sage: F*phi(X) == phi(f(X))
 210         True
 211
 212     """
 213     if is_MatrixSpace(M):
 214         basis_space = M
 215         basis = list(M.basis())
 216     else:
 217         basis_space = M[0].matrix_space()
 218         basis = M
 219
 220     (phi, phi_inv) = basis_repr(M)
 221
 222     # Get a basis for the image space. Since phi is an isometry,
 223     # it takes one basis to another.
 224     image_basis = [ phi(b) for b in basis ]
 225
 226     # Now construct the image space itself equipped with our custom basis.
 227     W = VectorSpace(basis_space.base_ring(), len(basis))
 228     W = W.span_of_basis(image_basis)
 229
 230     return matrix.column( W.coordinates(phi(L(b))) for b in basis )