From: Michael Orlitzky Date: Tue, 29 Oct 2019 01:57:11 +0000 (-0400) Subject: eja: drop custom gram_schmidt() routine that isn't noticeably faster. X-Git-Url: http://gitweb.michael.orlitzky.com/?p=sage.d.git;a=commitdiff_plain;h=1e9700cdd04434465ffcad148d078f7fa361e426 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 0f64141..c372e50 100644 --- 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 cf75e32..b6b0a03 100644 --- 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