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 _all2list(x):
   6     r"""
   7     Flatten a vector, matrix, or cartesian product of those things
   8     into a long list.
   9
  10     EXAMPLES::
  11
  12         sage: from mjo.eja.eja_utils import _all2list
  13         sage: V1 = VectorSpace(QQ,2)
  14         sage: V2 = MatrixSpace(QQ,2)
  15         sage: x1 = V1([1,1])
  16         sage: x2 = V1([1,-1])
  17         sage: y1 = V2.one()
  18         sage: y2 = V2([0,1,1,0])
  19         sage: _all2list((x1,y1))
  20         [1, 1, 1, 0, 0, 1]
  21         sage: _all2list((x2,y2))
  22         [1, -1, 0, 1, 1, 0]
  23         sage: M = cartesian_product([V1,V2])
  24         sage: _all2list(M((x1,y1)))
  25         [1, 1, 1, 0, 0, 1]
  26         sage: _all2list(M((x2,y2)))
  27         [1, -1, 0, 1, 1, 0]
  28
  29     """
  30     if hasattr(x, 'list'):
  31         # Easy case...
  32         return x.list()
  33     else:
  34         # But what if it's a tuple or something else? This has to
  35         # handle cartesian products of cartesian products, too; that's
  36         # why it's recursive.
  37         return sum( map(_all2list,x), [] )
  38
  39 def _mat2vec(m):
  40         return vector(m.base_ring(), m.list())
  41
  42 def _vec2mat(v):
  43         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
  44
  45 def gram_schmidt(v, inner_product=None):
  46     """
  47     Perform Gram-Schmidt on the list ``v`` which are assumed to be
  48     vectors over the same base ring. Returns a list of orthonormalized
  49     vectors over the smallest extention ring containing the necessary
  50     roots.
  51
  52     SETUP::
  53
  54         sage: from mjo.eja.eja_utils import gram_schmidt
  55
  56     EXAMPLES:
  57
  58     The usual inner-product and norm are default::
  59
  60         sage: v1 = vector(QQ,(1,2,3))
  61         sage: v2 = vector(QQ,(1,-1,6))
  62         sage: v3 = vector(QQ,(2,1,-1))
  63         sage: v = [v1,v2,v3]
  64         sage: u = gram_schmidt(v)
  65         sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
  66         True
  67         sage: bool(u[0].inner_product(u[1]) == 0)
  68         True
  69         sage: bool(u[0].inner_product(u[2]) == 0)
  70         True
  71         sage: bool(u[1].inner_product(u[2]) == 0)
  72         True
  73
  74
  75     But if you supply a custom inner product, the result is
  76     orthonormal with respect to that (and not the usual inner
  77     product)::
  78
  79         sage: v1 = vector(QQ,(1,2,3))
  80         sage: v2 = vector(QQ,(1,-1,6))
  81         sage: v3 = vector(QQ,(2,1,-1))
  82         sage: v = [v1,v2,v3]
  83         sage: B = matrix(QQ, [ [6, 4, 2],
  84         ....:                  [4, 5, 4],
  85         ....:                  [2, 4, 9] ])
  86         sage: ip = lambda x,y: (B*x).inner_product(y)
  87         sage: norm = lambda x: ip(x,x)
  88         sage: u = gram_schmidt(v,ip)
  89         sage: all( norm(u_i) == 1 for u_i in u )
  90         True
  91         sage: ip(u[0],u[1]).is_zero()
  92         True
  93         sage: ip(u[0],u[2]).is_zero()
  94         True
  95         sage: ip(u[1],u[2]).is_zero()
  96         True
  97
  98     This Gram-Schmidt routine can be used on matrices as well, so long
  99     as an appropriate inner-product is provided::
 100
 101         sage: E11 = matrix(QQ, [ [1,0],
 102         ....:                    [0,0] ])
 103         sage: E12 = matrix(QQ, [ [0,1],
 104         ....:                    [1,0] ])
 105         sage: E22 = matrix(QQ, [ [0,0],
 106         ....:                    [0,1] ])
 107         sage: I = matrix.identity(QQ,2)
 108         sage: trace_ip = lambda X,Y: (X*Y).trace()
 109         sage: gram_schmidt([E11,E12,I,E22], inner_product=trace_ip)
 110         [
 111         [1 0]  [          0 1/2*sqrt(2)]  [0 0]
 112         [0 0], [1/2*sqrt(2)           0], [0 1]
 113         ]
 114
 115     It even works on Cartesian product spaces whose factors are vector
 116     or matrix spaces::
 117
 118         sage: V1 = VectorSpace(AA,2)
 119         sage: V2 = MatrixSpace(AA,2)
 120         sage: M = cartesian_product([V1,V2])
 121         sage: x1 = V1([1,1])
 122         sage: x2 = V1([1,-1])
 123         sage: y1 = V2.one()
 124         sage: y2 = V2([0,1,1,0])
 125         sage: z1 = M((x1,y1))
 126         sage: z2 = M((x2,y2))
 127         sage: def ip(a,b):
 128         ....:     return a[0].inner_product(b[0]) + (a[1]*b[1]).trace()
 129         sage: U = gram_schmidt([z1,z2], inner_product=ip)
 130         sage: ip(U[0],U[1])
 131         0
 132         sage: ip(U[0],U[0])
 133         1
 134         sage: ip(U[1],U[1])
 135         1
 136
 137     TESTS:
 138
 139     Ensure that zero vectors don't get in the way::
 140
 141         sage: v1 = vector(QQ,(1,2,3))
 142         sage: v2 = vector(QQ,(1,-1,6))
 143         sage: v3 = vector(QQ,(0,0,0))
 144         sage: v = [v1,v2,v3]
 145         sage: len(gram_schmidt(v)) == 2
 146         True
 147     """
 148     if inner_product is None:
 149         inner_product = lambda x,y: x.inner_product(y)
 150     norm = lambda x: inner_product(x,x).sqrt()
 151
 152     v = list(v) # make a copy, don't clobber the input
 153
 154     # Drop all zero vectors before we start.
 155     v = [ v_i for v_i in v if not v_i.is_zero() ]
 156
 157     if len(v) == 0:
 158         # cool
 159         return v
 160
 161     R = v[0].base_ring()
 162
 163     # Define a scaling operation that can be used on tuples.
 164     # Oh and our "zero" needs to belong to the right space.
 165     scale = lambda x,alpha: x*alpha
 166     zero = v[0].parent().zero()
 167     if hasattr(v[0], 'cartesian_factors'):
 168         P = v[0].parent()
 169         scale = lambda x,alpha: P(tuple( x_i*alpha
 170                                          for x_i in x.cartesian_factors() ))
 171
 172
 173     def proj(x,y):
 174         return scale(x, (inner_product(x,y)/inner_product(x,x)))
 175
 176     # First orthogonalize...
 177     for i in range(1,len(v)):
 178         # Earlier vectors can be made into zero so we have to ignore them.
 179         v[i] -= sum( (proj(v[j],v[i])
 180                       for j in range(i)
 181                       if not v[j].is_zero() ),
 182                      zero )
 183
 184     # And now drop all zero vectors again if they were "orthogonalized out."
 185     v = [ v_i for v_i in v if not v_i.is_zero() ]
 186
 187     # Just normalize. If the algebra is missing the roots, we can't add
 188     # them here because then our subalgebra would have a bigger field
 189     # than the superalgebra.
 190     for i in range(len(v)):
 191         v[i] = scale(v[i], ~norm(v[i]))
 192
 193     return v