mjo/eja/eja_utils.py

   1 from sage.structure.element import is_Matrix
   2
   3 def _scale(x, alpha):
   4     r"""
   5     Scale the vector, matrix, or cartesian-product-of-those-things
   6     ``x`` by ``alpha``.
   7
   8     This works around the inability to scale certain elements of
   9     Cartesian product spaces, as reported in
  10
  11       https://trac.sagemath.org/ticket/31435
  12
  13     ..WARNING:
  14
  15         This will do the wrong thing if you feed it a tuple or list.
  16
  17     SETUP::
  18
  19         sage: from mjo.eja.eja_utils import _scale
  20
  21     EXAMPLES::
  22
  23         sage: v = vector(QQ, (1,2,3))
  24         sage: _scale(v,2)
  25         (2, 4, 6)
  26         sage: m = matrix(QQ, [[1,2],[3,4]])
  27         sage: M = cartesian_product([m.parent(), m.parent()])
  28         sage: _scale(M((m,m)), 2)
  29         ([2 4]
  30         [6 8], [2 4]
  31         [6 8])
  32
  33     """
  34     if hasattr(x, 'cartesian_factors'):
  35         P = x.parent()
  36         return P(tuple( _scale(x_i, alpha)
  37                         for x_i in x.cartesian_factors() ))
  38     else:
  39         return x*alpha
  40
  41
  42 def _all2list(x):
  43     r"""
  44     Flatten a vector, matrix, or cartesian product of those things
  45     into a long list.
  46
  47     If the entries of the matrix themselves belong to a real vector
  48     space (such as the complex numbers which can be thought of as
  49     pairs of real numbers), they will also be expanded in vector form
  50     and flattened into the list.
  51
  52     SETUP::
  53
  54         sage: from mjo.eja.eja_utils import _all2list
  55         sage: from mjo.hurwitz import (QuaternionMatrixAlgebra,
  56         ....:                          Octonions,
  57         ....:                          OctonionMatrixAlgebra)
  58
  59     EXAMPLES::
  60
  61         sage: _all2list([[1]])
  62         [1]
  63
  64     ::
  65
  66         sage: V1 = VectorSpace(QQ,2)
  67         sage: V2 = MatrixSpace(QQ,2)
  68         sage: x1 = V1([1,1])
  69         sage: x2 = V1([1,-1])
  70         sage: y1 = V2.one()
  71         sage: y2 = V2([0,1,1,0])
  72         sage: _all2list((x1,y1))
  73         [1, 1, 1, 0, 0, 1]
  74         sage: _all2list((x2,y2))
  75         [1, -1, 0, 1, 1, 0]
  76         sage: M = cartesian_product([V1,V2])
  77         sage: _all2list(M((x1,y1)))
  78         [1, 1, 1, 0, 0, 1]
  79         sage: _all2list(M((x2,y2)))
  80         [1, -1, 0, 1, 1, 0]
  81
  82     ::
  83
  84         sage: _all2list(Octonions().one())
  85         [1, 0, 0, 0, 0, 0, 0, 0]
  86         sage: _all2list(OctonionMatrixAlgebra(1).one())
  87         [1, 0, 0, 0, 0, 0, 0, 0]
  88
  89     ::
  90
  91         sage: _all2list(QuaternionAlgebra(QQ, -1, -1).one())
  92         [1, 0, 0, 0]
  93         sage: _all2list(QuaternionMatrixAlgebra(1).one())
  94         [1, 0, 0, 0]
  95
  96     ::
  97
  98         sage: V1 = VectorSpace(QQ,2)
  99         sage: V2 = OctonionMatrixAlgebra(1,field=QQ)
 100         sage: C = cartesian_product([V1,V2])
 101         sage: x1 = V1([3,4])
 102         sage: y1 = V2.one()
 103         sage: _all2list(C( (x1,y1) ))
 104         [3, 4, 1, 0, 0, 0, 0, 0, 0, 0]
 105
 106     """
 107     if hasattr(x, 'to_vector'):
 108         # This works on matrices of e.g. octonions directly, without
 109         # first needing to convert them to a list of octonions and
 110         # then recursing down into the list. It also avoids the wonky
 111         # list(x) when x is an element of a CFM. I don't know what it
 112         # returns but it aint the coordinates. We don't recurse
 113         # because vectors can only contain ring elements as entries.
 114         return x.to_vector().list()
 115
 116     if is_Matrix(x):
 117         # This sucks, but for performance reasons we don't want to
 118         # call _all2list recursively on the contents of a matrix
 119         # when we don't have to (they only contain ring elements
 120         # as entries)
 121         return x.list()
 122
 123     try:
 124         xl = list(x)
 125     except TypeError: # x is not iterable
 126         return [x]
 127
 128     if xl == [x]:
 129         # Avoid the retardation of list(QQ(1)) == [1].
 130         return [x]
 131
 132     return sum( map(_all2list, xl) , [])
 133
 134
 135 def gram_schmidt(v, inner_product=None):
 136     """
 137     Perform Gram-Schmidt on the list ``v`` which are assumed to be
 138     vectors over the same base ring. Returns a list of orthonormalized
 139     vectors over the same base ring, which means that your base ring
 140     needs to contain the appropriate roots.
 141
 142     SETUP::
 143
 144         sage: from mjo.eja.eja_utils import gram_schmidt
 145
 146     EXAMPLES:
 147
 148     If you start with an orthonormal set, you get it back. We can use
 149     the rationals here because we don't need any square roots::
 150
 151         sage: v1 = vector(QQ, (1,0,0))
 152         sage: v2 = vector(QQ, (0,1,0))
 153         sage: v3 = vector(QQ, (0,0,1))
 154         sage: v = [v1,v2,v3]
 155         sage: gram_schmidt(v) == v
 156         True
 157
 158     The usual inner-product and norm are default::
 159
 160         sage: v1 = vector(AA,(1,2,3))
 161         sage: v2 = vector(AA,(1,-1,6))
 162         sage: v3 = vector(AA,(2,1,-1))
 163         sage: v = [v1,v2,v3]
 164         sage: u = gram_schmidt(v)
 165         sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
 166         True
 167         sage: bool(u[0].inner_product(u[1]) == 0)
 168         True
 169         sage: bool(u[0].inner_product(u[2]) == 0)
 170         True
 171         sage: bool(u[1].inner_product(u[2]) == 0)
 172         True
 173
 174
 175     But if you supply a custom inner product, the result is
 176     orthonormal with respect to that (and not the usual inner
 177     product)::
 178
 179         sage: v1 = vector(AA,(1,2,3))
 180         sage: v2 = vector(AA,(1,-1,6))
 181         sage: v3 = vector(AA,(2,1,-1))
 182         sage: v = [v1,v2,v3]
 183         sage: B = matrix(AA, [ [6, 4, 2],
 184         ....:                  [4, 5, 4],
 185         ....:                  [2, 4, 9] ])
 186         sage: ip = lambda x,y: (B*x).inner_product(y)
 187         sage: norm = lambda x: ip(x,x)
 188         sage: u = gram_schmidt(v,ip)
 189         sage: all( norm(u_i) == 1 for u_i in u )
 190         True
 191         sage: ip(u[0],u[1]).is_zero()
 192         True
 193         sage: ip(u[0],u[2]).is_zero()
 194         True
 195         sage: ip(u[1],u[2]).is_zero()
 196         True
 197
 198     This Gram-Schmidt routine can be used on matrices as well, so long
 199     as an appropriate inner-product is provided::
 200
 201         sage: E11 = matrix(AA, [ [1,0],
 202         ....:                    [0,0] ])
 203         sage: E12 = matrix(AA, [ [0,1],
 204         ....:                    [1,0] ])
 205         sage: E22 = matrix(AA, [ [0,0],
 206         ....:                    [0,1] ])
 207         sage: I = matrix.identity(AA,2)
 208         sage: trace_ip = lambda X,Y: (X*Y).trace()
 209         sage: gram_schmidt([E11,E12,I,E22], inner_product=trace_ip)
 210         [
 211         [1 0]  [                  0 0.7071067811865475?]  [0 0]
 212         [0 0], [0.7071067811865475?                   0], [0 1]
 213         ]
 214
 215     It even works on Cartesian product spaces whose factors are vector
 216     or matrix spaces::
 217
 218         sage: V1 = VectorSpace(AA,2)
 219         sage: V2 = MatrixSpace(AA,2)
 220         sage: M = cartesian_product([V1,V2])
 221         sage: x1 = V1([1,1])
 222         sage: x2 = V1([1,-1])
 223         sage: y1 = V2.one()
 224         sage: y2 = V2([0,1,1,0])
 225         sage: z1 = M((x1,y1))
 226         sage: z2 = M((x2,y2))
 227         sage: def ip(a,b):
 228         ....:     return a[0].inner_product(b[0]) + (a[1]*b[1]).trace()
 229         sage: U = gram_schmidt([z1,z2], inner_product=ip)
 230         sage: ip(U[0],U[1])
 231         0
 232         sage: ip(U[0],U[0])
 233         1
 234         sage: ip(U[1],U[1])
 235         1
 236
 237     TESTS:
 238
 239     Ensure that zero vectors don't get in the way::
 240
 241         sage: v1 = vector(AA,(1,2,3))
 242         sage: v2 = vector(AA,(1,-1,6))
 243         sage: v3 = vector(AA,(0,0,0))
 244         sage: v = [v1,v2,v3]
 245         sage: len(gram_schmidt(v)) == 2
 246         True
 247     """
 248     if len(v) == 0:
 249         # cool
 250         return v
 251
 252     V = v[0].parent()
 253
 254     if inner_product is None:
 255         inner_product = lambda x,y: x.inner_product(y)
 256
 257     def norm(x):
 258         # Don't expand the given field; the inner-product's codomain
 259         # is already correct. For example QQ(2).sqrt() returns sqrt(2)
 260         # in SR, and that will give you weird errors about symbolics
 261         # when what's really going wrong is that you're trying to
 262         # orthonormalize in QQ.
 263         return V.base_ring()(inner_product(x,x).sqrt())
 264
 265     sc = lambda x,a: a*x
 266     if hasattr(V, 'cartesian_factors'):
 267         # Only use the slow implementation if necessary.
 268         sc = _scale
 269
 270     def proj(x,y):
 271         # project y onto the span of {x}
 272         return sc(x, (inner_product(x,y)/inner_product(x,x)))
 273
 274     def normalize(x):
 275         return sc(x, ~norm(x))
 276
 277     v_out = [] # make a copy, don't clobber the input
 278
 279     for (i, v_i) in enumerate(v):
 280         ortho_v_i = v_i - V.sum( proj(v_out[j],v_i) for j in range(i) )
 281         if not ortho_v_i.is_zero():
 282             v_out.append(normalize(ortho_v_i))
 283
 284     return v_out