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 from sage.rings.number_field.number_field import NumberField
   5 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
   6 from sage.rings.real_lazy import RLF
   7
   8 def _mat2vec(m):
   9         return vector(m.base_ring(), m.list())
  10
  11 def _vec2mat(v):
  12         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
  13
  14 def gram_schmidt(v, inner_product=None):
  15     """
  16     Perform Gram-Schmidt on the list ``v`` which are assumed to be
  17     vectors over the same base ring. Returns a list of orthonormalized
  18     vectors over the smallest extention ring containing the necessary
  19     roots.
  20
  21     SETUP::
  22
  23         sage: from mjo.eja.eja_utils import gram_schmidt
  24
  25     EXAMPLES:
  26
  27     The usual inner-product and norm are default::
  28
  29         sage: v1 = vector(QQ,(1,2,3))
  30         sage: v2 = vector(QQ,(1,-1,6))
  31         sage: v3 = vector(QQ,(2,1,-1))
  32         sage: v = [v1,v2,v3]
  33         sage: u = gram_schmidt(v)
  34         sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
  35         True
  36         sage: bool(u[0].inner_product(u[1]) == 0)
  37         True
  38         sage: bool(u[0].inner_product(u[2]) == 0)
  39         True
  40         sage: bool(u[1].inner_product(u[2]) == 0)
  41         True
  42
  43
  44     But if you supply a custom inner product, the result is
  45     orthonormal with respect to that (and not the usual inner
  46     product)::
  47
  48         sage: v1 = vector(QQ,(1,2,3))
  49         sage: v2 = vector(QQ,(1,-1,6))
  50         sage: v3 = vector(QQ,(2,1,-1))
  51         sage: v = [v1,v2,v3]
  52         sage: B = matrix(QQ, [ [6, 4, 2],
  53         ....:                  [4, 5, 4],
  54         ....:                  [2, 4, 9] ])
  55         sage: ip = lambda x,y: (B*x).inner_product(y)
  56         sage: norm = lambda x: ip(x,x)
  57         sage: u = gram_schmidt(v,ip)
  58         sage: all( norm(u_i) == 1 for u_i in u )
  59         True
  60         sage: ip(u[0],u[1]).is_zero()
  61         True
  62         sage: ip(u[0],u[2]).is_zero()
  63         True
  64         sage: ip(u[1],u[2]).is_zero()
  65         True
  66
  67     TESTS:
  68
  69     Ensure that zero vectors don't get in the way::
  70
  71         sage: v1 = vector(QQ,(1,2,3))
  72         sage: v2 = vector(QQ,(1,-1,6))
  73         sage: v3 = vector(QQ,(0,0,0))
  74         sage: v = [v1,v2,v3]
  75         sage: len(gram_schmidt(v)) == 2
  76         True
  77
  78     """
  79     if inner_product is None:
  80         inner_product = lambda x,y: x.inner_product(y)
  81     norm = lambda x: inner_product(x,x).sqrt()
  82
  83     def proj(x,y):
  84         return (inner_product(x,y)/inner_product(x,x))*x
  85
  86     v = list(v) # make a copy, don't clobber the input
  87
  88     # Drop all zero vectors before we start.
  89     v = [ v_i for v_i in v if not v_i.is_zero() ]
  90
  91     if len(v) == 0:
  92         # cool
  93         return v
  94
  95     R = v[0].base_ring()
  96
  97     # First orthogonalize...
  98     for i in range(1,len(v)):
  99         # Earlier vectors can be made into zero so we have to ignore them.
 100         v[i] -= sum( proj(v[j],v[i]) for j in range(i) if not v[j].is_zero() )
 101
 102     # And now drop all zero vectors again if they were "orthogonalized out."
 103     v = [ v_i for v_i in v if not v_i.is_zero() ]
 104
 105     # Just normalize. If the algebra is missing the roots, we can't add
 106     # them here because then our subalgebra would have a bigger field
 107     # than the superalgebra.
 108     for i in range(len(v)):
 109         v[i] = v[i] / norm(v[i])
 110
 111     return v