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