From f72c84ce3d46f2611a65417c72e9017754ec156f Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Thu, 29 Aug 2019 08:56:56 -0400 Subject: [PATCH] eja: add a WIP gram-schmidt for EJA elements. This doesn't really work right now because we need a whole bunch of algebraic numbers that we don't know a priori. I might need to suck it up and just use AA instead of quadratic number fields. --- mjo/eja/eja_element.py | 4 +- mjo/eja/eja_subalgebra.py | 44 +++++++++++-------- mjo/eja/eja_utils.py | 91 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 120 insertions(+), 19 deletions(-) diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py index 0705018..c2f9fe6 100644 --- a/mjo/eja/eja_element.py +++ b/mjo/eja/eja_element.py @@ -1011,7 +1011,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): - def subalgebra_generated_by(self): + def subalgebra_generated_by(self, orthonormalize_basis=False): """ Return the associative subalgebra of the parent EJA generated by this element. @@ -1050,7 +1050,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 0 """ - return FiniteDimensionalEuclideanJordanElementSubalgebra(self) + return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis) def subalgebra_idempotent(self): diff --git a/mjo/eja/eja_subalgebra.py b/mjo/eja/eja_subalgebra.py index 8f6e56b..fee718b 100644 --- a/mjo/eja/eja_subalgebra.py +++ b/mjo/eja/eja_subalgebra.py @@ -2,7 +2,7 @@ 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): """ @@ -99,7 +99,7 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide 1 """ - def __init__(self, elt): + def __init__(self, elt, orthonormalize_basis): self._superalgebra = elt.parent() category = self._superalgebra.category().Associative() V = self._superalgebra.vector_space() @@ -135,26 +135,36 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide natural_basis=natural_basis) - # First compute the vector subspace spanned by the powers of - # the given element. + # This list is guaranteed to contain all independent powers, + # 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 ] - # Figure out which powers form a linearly-independent set. - ind_rows = matrix(field, power_vectors).pivot_rows() + 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 ] - # Pick those out of the list of all powers. - superalgebra_basis = tuple(map(powers.__getitem__, ind_rows)) + # Figure out which powers form a linearly-independent set. + ind_rows = matrix(field, power_vectors).pivot_rows() - # If our superalgebra is a subalgebra of something else, then - # these vectors won't have the right coordinates for - # V.span_of_basis() unless we use V.from_vector() on them. - basis_vectors = map(power_vectors.__getitem__, ind_rows) - W = V.span_of_basis( V.from_vector(v) for v in basis_vectors ) + # Pick those out of the list of all powers. + superalgebra_basis = tuple(map(powers.__getitem__, ind_rows)) - # Now figure out the entries of the right-multiplication - # matrix for the successive basis elements b0, b1,... of - # that subspace. + # If our superalgebra is a subalgebra of something else, then + # these vectors won't have the right coordinates for + # V.span_of_basis() unless we use V.from_vector() on them. + basis_vectors = map(power_vectors.__getitem__, ind_rows) + 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 ] + + W = V.span_of_basis( V.from_vector(v) for v in basis_vectors ) n = len(superalgebra_basis) mult_table = [[W.zero() for i in range(n)] for j in range(n)] for i in range(n): diff --git a/mjo/eja/eja_utils.py b/mjo/eja/eja_utils.py index b6b0a03..8f2d8f3 100644 --- a/mjo/eja/eja_utils.py +++ b/mjo/eja/eja_utils.py @@ -1,4 +1,95 @@ 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: [ u_i.inner_product(u_i).sqrt() == 1 for u_i in u ] + True + sage: u[0].inner_product(u[1]) == 0 + True + sage: u[0].inner_product(u[2]) == 0 + True + sage: 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() ] + + # Now pretend to normalize, building a new ring R that contains + # all of the necessary square roots. + norms_squared = [0]*len(v) + + for i in xrange(len(v)): + norms_squared[i] = v[i].inner_product(v[i]) + ns = [norms_squared[i].numerator(), norms_squared[i].denominator()] + + # Do the numerator and denominator separately so that we + # adjoin e.g. sqrt(2) and sqrt(3) instead of sqrt(2/3). + for j in xrange(len(ns)): + PR = PolynomialRing(R, 'z') + z = PR.gen() + p = z**2 - ns[j] + if p.is_irreducible(): + R = NumberField(p, + 'sqrt' + str(ns[j]), + embedding=RLF(ns[j]).sqrt()) + + # When we're done, we have to change every element's ring to the + # extension that we wound up with, and then normalize it (which + # should work, since "R" contains its norm now). + for i in xrange(len(v)): + v[i] = v[i].change_ring(R) / R(norms_squared[i]).sqrt() + + return v -- 2.44.2