eja: drop custom gram_schmidt() routine that isn't noticeably faster.

Using Gram-Schmidt in SageMath is a bit of a pain, since you have to
put your vector space elements into a matrix, which means you have to
convert them to vectors first. And then you wind up transposing the
matrix twice, and then you have to convert everything back to vectors
and then back to algebra elements, dropping the zero vectors/elements
in the process. Sounds slow, right?

Well, testing shows that jumping through all those hoops and using the
default implementation isn't much slower than the custom routine I
wrote. So, for now, let's delete it.

Beware that this breaks the construction of subalgebras when the field
is not exact because SageMath won't orthonormalize over an inexact
field.

diff --git a/mjo/eja/eja_subalgebra.py b/mjo/eja/eja_subalgebra.py
index 0f641416366ef7ee72d4703e070c69e2d60e30a9..c372e50072c73a50281dd45d07eec62a62848277 100644 (file)
--- a/mjo/eja/eja_subalgebra.py
+++ b/mjo/eja/eja_subalgebra.py
@@ -2,7 +2,6 @@ from sage.matrix.constructor import matrix
 
 from mjo.eja.eja_algebra import FiniteDimensionalEuclideanJordanAlgebra
 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
-from mjo.eja.eja_utils import gram_schmidt
 
 class FiniteDimensionalEuclideanJordanElementSubalgebraElement(FiniteDimensionalEuclideanJordanAlgebraElement):
     """
@@ -132,16 +131,17 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         # because it's the maximal set of powers that could possibly
         # be independent (by a dimension argument).
         powers = [ elt**k for k in range(V.dimension()) ]
+        power_vectors = [ p.to_vector() for p in powers ]
+        P = matrix(field, power_vectors)
 
         if orthonormalize_basis == False:
             # In this case, we just need to figure out which elements
             # of the "powers" list are redundant... First compute the
             # vector subspace spanned by the powers of the given
             # element.
-            power_vectors = [ p.to_vector() for p in powers ]
 
             # Figure out which powers form a linearly-independent set.
-            ind_rows = matrix(field, power_vectors).pivot_rows()
+            ind_rows = P.pivot_rows()
 
             # Pick those out of the list of all powers.
             superalgebra_basis = tuple(map(powers.__getitem__, ind_rows))
@@ -153,9 +153,16 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         else:
             # If we're going to orthonormalize the basis anyway, we
             # might as well just do Gram-Schmidt on the whole list of
-            # powers. The redundant ones will get zero'd out.
-            superalgebra_basis = gram_schmidt(powers)
-            basis_vectors = [ b.to_vector() for b in superalgebra_basis ]
+            # powers. The redundant ones will get zero'd out. If this
+            # looks like a roundabout way to orthonormalize, it is.
+            # But converting everything from algebra elements to vectors
+            # to matrices and then back again turns out to be about
+            # as fast as reimplementing our own Gram-Schmidt that
+            # works in an EJA.
+            G,_ = P.gram_schmidt(orthonormal=True)
+            basis_vectors = [ g for g in G.rows() if not g.is_zero() ]
+            superalgebra_basis = [ self._superalgebra.from_vector(b)
+                                   for b in basis_vectors ]
 
         W = V.span_of_basis( V.from_vector(v) for v in basis_vectors )
         n = len(superalgebra_basis)

diff --git a/mjo/eja/eja_utils.py b/mjo/eja/eja_utils.py
index cf75e325697dcefb3bf682b855f8d83e3e4f89e2..b6b0a0327c19842ce654023b00bf6e0b9af539a1 100644 (file)
--- a/mjo/eja/eja_utils.py
+++ b/mjo/eja/eja_utils.py
@@ -1,76 +1,4 @@
 from sage.modules.free_module_element import vector
-from sage.rings.number_field.number_field import NumberField
-from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-from sage.rings.real_lazy import RLF
 
 def _mat2vec(m):
         return vector(m.base_ring(), m.list())
-
-def gram_schmidt(v):
-    """
-    Perform Gram-Schmidt on the list ``v`` which are assumed to be
-    vectors over the same base ring. Returns a list of orthonormalized
-    vectors over the smallest extention ring containing the necessary
-    roots.
-
-    SETUP::
-
-        sage: from mjo.eja.eja_utils import gram_schmidt
-
-    EXAMPLES::
-
-        sage: v1 = vector(QQ,(1,2,3))
-        sage: v2 = vector(QQ,(1,-1,6))
-        sage: v3 = vector(QQ,(2,1,-1))
-        sage: v = [v1,v2,v3]
-        sage: u = gram_schmidt(v)
-        sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
-        True
-        sage: bool(u[0].inner_product(u[1]) == 0)
-        True
-        sage: bool(u[0].inner_product(u[2]) == 0)
-        True
-        sage: bool(u[1].inner_product(u[2]) == 0)
-        True
-
-    TESTS:
-
-    Ensure that zero vectors don't get in the way::
-
-        sage: v1 = vector(QQ,(1,2,3))
-        sage: v2 = vector(QQ,(1,-1,6))
-        sage: v3 = vector(QQ,(0,0,0))
-        sage: v = [v1,v2,v3]
-        sage: len(gram_schmidt(v)) == 2
-        True
-
-    """
-    def proj(x,y):
-        return (y.inner_product(x)/x.inner_product(x))*x
-
-    v = list(v) # make a copy, don't clobber the input
-
-    # Drop all zero vectors before we start.
-    v = [ v_i for v_i in v if not v_i.is_zero() ]
-
-    if len(v) == 0:
-        # cool
-        return v
-
-    R = v[0].base_ring()
-
-    # First orthogonalize...
-    for i in xrange(1,len(v)):
-        # Earlier vectors can be made into zero so we have to ignore them.
-        v[i] -= sum( proj(v[j],v[i]) for j in range(i) if not v[j].is_zero() )
-
-    # And now drop all zero vectors again if they were "orthogonalized out."
-    v = [ v_i for v_i in v if not v_i.is_zero() ]
-
-    # Just normalize. If the algebra is missing the roots, we can't add
-    # them here because then our subalgebra would have a bigger field
-    # than the superalgebra.
-    for i in xrange(len(v)):
-        v[i] = v[i] / v[i].norm()
-
-    return v