mjo/eja/eja_utils.py

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