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.octonions import Octonions, 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: V1 = VectorSpace(QQ,2)
  92         sage: V2 = OctonionMatrixAlgebra(1,field=QQ)
  93         sage: C = cartesian_product([V1,V2])
  94         sage: x1 = V1([3,4])
  95         sage: y1 = V2.one()
  96         sage: _all2list(C( (x1,y1) ))
  97         [3, 4, 1, 0, 0, 0, 0, 0, 0, 0]
  98
  99     """
 100     if hasattr(x, 'to_vector'):
 101         # This works on matrices of e.g. octonions directly, without
 102         # first needing to convert them to a list of octonions and
 103         # then recursing down into the list. It also avoids the wonky
 104         # list(x) when x is an element of a CFM. I don't know what it
 105         # returns but it aint the coordinates. This will fall through
 106         # to the iterable case the next time around.
 107         return _all2list(x.to_vector())
 108
 109     try:
 110         xl = list(x)
 111     except TypeError: # x is not iterable
 112         return [x]
 113
 114     if xl == [x]:
 115         # Avoid the retardation of list(QQ(1)) == [1].
 116         return [x]
 117
 118     return sum(list( map(_all2list, xl) ), [])
 119
 120
 121
 122 def _mat2vec(m):
 123         return vector(m.base_ring(), m.list())
 124
 125 def _vec2mat(v):
 126         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
 127
 128 def gram_schmidt(v, inner_product=None):
 129     """
 130     Perform Gram-Schmidt on the list ``v`` which are assumed to be
 131     vectors over the same base ring. Returns a list of orthonormalized
 132     vectors over the smallest extention ring containing the necessary
 133     roots.
 134
 135     SETUP::
 136
 137         sage: from mjo.eja.eja_utils import gram_schmidt
 138
 139     EXAMPLES:
 140
 141     The usual inner-product and norm are default::
 142
 143         sage: v1 = vector(QQ,(1,2,3))
 144         sage: v2 = vector(QQ,(1,-1,6))
 145         sage: v3 = vector(QQ,(2,1,-1))
 146         sage: v = [v1,v2,v3]
 147         sage: u = gram_schmidt(v)
 148         sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
 149         True
 150         sage: bool(u[0].inner_product(u[1]) == 0)
 151         True
 152         sage: bool(u[0].inner_product(u[2]) == 0)
 153         True
 154         sage: bool(u[1].inner_product(u[2]) == 0)
 155         True
 156
 157
 158     But if you supply a custom inner product, the result is
 159     orthonormal with respect to that (and not the usual inner
 160     product)::
 161
 162         sage: v1 = vector(QQ,(1,2,3))
 163         sage: v2 = vector(QQ,(1,-1,6))
 164         sage: v3 = vector(QQ,(2,1,-1))
 165         sage: v = [v1,v2,v3]
 166         sage: B = matrix(QQ, [ [6, 4, 2],
 167         ....:                  [4, 5, 4],
 168         ....:                  [2, 4, 9] ])
 169         sage: ip = lambda x,y: (B*x).inner_product(y)
 170         sage: norm = lambda x: ip(x,x)
 171         sage: u = gram_schmidt(v,ip)
 172         sage: all( norm(u_i) == 1 for u_i in u )
 173         True
 174         sage: ip(u[0],u[1]).is_zero()
 175         True
 176         sage: ip(u[0],u[2]).is_zero()
 177         True
 178         sage: ip(u[1],u[2]).is_zero()
 179         True
 180
 181     This Gram-Schmidt routine can be used on matrices as well, so long
 182     as an appropriate inner-product is provided::
 183
 184         sage: E11 = matrix(QQ, [ [1,0],
 185         ....:                    [0,0] ])
 186         sage: E12 = matrix(QQ, [ [0,1],
 187         ....:                    [1,0] ])
 188         sage: E22 = matrix(QQ, [ [0,0],
 189         ....:                    [0,1] ])
 190         sage: I = matrix.identity(QQ,2)
 191         sage: trace_ip = lambda X,Y: (X*Y).trace()
 192         sage: gram_schmidt([E11,E12,I,E22], inner_product=trace_ip)
 193         [
 194         [1 0]  [          0 1/2*sqrt(2)]  [0 0]
 195         [0 0], [1/2*sqrt(2)           0], [0 1]
 196         ]
 197
 198     It even works on Cartesian product spaces whose factors are vector
 199     or matrix spaces::
 200
 201         sage: V1 = VectorSpace(AA,2)
 202         sage: V2 = MatrixSpace(AA,2)
 203         sage: M = cartesian_product([V1,V2])
 204         sage: x1 = V1([1,1])
 205         sage: x2 = V1([1,-1])
 206         sage: y1 = V2.one()
 207         sage: y2 = V2([0,1,1,0])
 208         sage: z1 = M((x1,y1))
 209         sage: z2 = M((x2,y2))
 210         sage: def ip(a,b):
 211         ....:     return a[0].inner_product(b[0]) + (a[1]*b[1]).trace()
 212         sage: U = gram_schmidt([z1,z2], inner_product=ip)
 213         sage: ip(U[0],U[1])
 214         0
 215         sage: ip(U[0],U[0])
 216         1
 217         sage: ip(U[1],U[1])
 218         1
 219
 220     TESTS:
 221
 222     Ensure that zero vectors don't get in the way::
 223
 224         sage: v1 = vector(QQ,(1,2,3))
 225         sage: v2 = vector(QQ,(1,-1,6))
 226         sage: v3 = vector(QQ,(0,0,0))
 227         sage: v = [v1,v2,v3]
 228         sage: len(gram_schmidt(v)) == 2
 229         True
 230     """
 231     if inner_product is None:
 232         inner_product = lambda x,y: x.inner_product(y)
 233     def norm(x):
 234         ip = inner_product(x,x)
 235         # Don't expand the given field; the inner-product's codomain
 236         # is already correct. For example QQ(2).sqrt() returns sqrt(2)
 237         # in SR, and that will give you weird errors about symbolics
 238         # when what's really going wrong is that you're trying to
 239         # orthonormalize in QQ.
 240         return ip.parent()(ip.sqrt())
 241
 242     v = list(v) # make a copy, don't clobber the input
 243
 244     # Drop all zero vectors before we start.
 245     v = [ v_i for v_i in v if not v_i.is_zero() ]
 246
 247     if len(v) == 0:
 248         # cool
 249         return v
 250
 251     R = v[0].base_ring()
 252
 253     # Our "zero" needs to belong to the right space for sum() to work.
 254     zero = v[0].parent().zero()
 255
 256     sc = lambda x,a: a*x
 257     if hasattr(v[0], 'cartesian_factors'):
 258         # Only use the slow implementation if necessary.
 259         sc = _scale
 260
 261     def proj(x,y):
 262         return sc(x, (inner_product(x,y)/inner_product(x,x)))
 263
 264     # First orthogonalize...
 265     for i in range(1,len(v)):
 266         # Earlier vectors can be made into zero so we have to ignore them.
 267         v[i] -= sum( (proj(v[j],v[i])
 268                       for j in range(i)
 269                       if not v[j].is_zero() ),
 270                      zero )
 271
 272     # And now drop all zero vectors again if they were "orthogonalized out."
 273     v = [ v_i for v_i in v if not v_i.is_zero() ]
 274
 275     # Just normalize. If the algebra is missing the roots, we can't add
 276     # them here because then our subalgebra would have a bigger field
 277     # than the superalgebra.
 278     for i in range(len(v)):
 279         v[i] = sc(v[i], ~norm(v[i]))
 280
 281     return v