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