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