From: Michael Orlitzky Date: Wed, 21 Aug 2019 20:10:14 +0000 (-0400) Subject: eja: refactor some of the basis and inner-product stuff. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=11e681d6320f0b7ddbb834931845b6f4a745da93;p=sage.d.git eja: refactor some of the basis and inner-product stuff. This is movement towards eventually cheating on the charpoly coefficients, which we should be able to compute in the "nice" basis and then scale to the normalized one. The coefficients are polynomials in "the coordinates of x", and those coordinates change only by a scalar multiple when we normalize the basis. --- diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index a7e5618..d0e7b07 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -416,9 +416,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): True """ - if (not x in self) or (not y in self): - raise TypeError("arguments must live in this algebra") - return x.trace_inner_product(y) + X = x.natural_representation() + Y = y.natural_representation() + return self.__class__.natural_inner_product(X,Y) def is_trivial(self): @@ -537,6 +537,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return self._natural_basis[0].matrix_space() + @staticmethod + def natural_inner_product(X,Y): + """ + Compute the inner product of two naturally-represented elements. + + For example in the real symmetric matrix EJA, this will compute + the trace inner-product of two n-by-n symmetric matrices. The + default should work for the real cartesian product EJA, the + Jordan spin EJA, and the real symmetric matrices. The others + will have to be overridden. + """ + return (X.conjugate_transpose()*Y).trace() + + @cached_method def one(self): """ @@ -748,7 +762,30 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra): return fdeja.__init__(field, mult_table, rank=n, **kwargs) def inner_product(self, x, y): - return _usual_ip(x,y) + """ + Faster to reimplement than to use natural representations. + + SETUP:: + + sage: from mjo.eja.eja_algebra import RealCartesianProductEJA + + TESTS: + + Ensure that this is the usual inner product for the algebras + over `R^n`:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: J = RealCartesianProductEJA(n) + sage: x = J.random_element() + sage: y = J.random_element() + sage: X = x.natural_representation() + sage: Y = y.natural_representation() + sage: x.inner_product(y) == J.__class__.natural_inner_product(X,Y) + True + + """ + return x.to_vector().inner_product(y.to_vector()) def random_eja(): @@ -802,7 +839,7 @@ def random_eja(): -def _real_symmetric_basis(n, field, normalize): +def _real_symmetric_basis(n, field): """ Return a basis for the space of real symmetric n-by-n matrices. @@ -814,7 +851,7 @@ def _real_symmetric_basis(n, field, normalize): sage: set_random_seed() sage: n = ZZ.random_element(1,5) - sage: B = _real_symmetric_basis(n, QQbar, False) + sage: B = _real_symmetric_basis(n, QQ) sage: all( M.is_symmetric() for M in B) True @@ -829,13 +866,11 @@ def _real_symmetric_basis(n, field, normalize): Sij = Eij else: Sij = Eij + Eij.transpose() - if normalize: - Sij = Sij / _real_symmetric_matrix_ip(Sij,Sij).sqrt() S.append(Sij) return tuple(S) -def _complex_hermitian_basis(n, field, normalize): +def _complex_hermitian_basis(n, field): """ Returns a basis for the space of complex Hermitian n-by-n matrices. @@ -854,7 +889,7 @@ def _complex_hermitian_basis(n, field, normalize): sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: field = QuadraticField(2, 'sqrt2') - sage: B = _complex_hermitian_basis(n, field, False) + sage: B = _complex_hermitian_basis(n, field) sage: all( M.is_symmetric() for M in B) True @@ -877,8 +912,7 @@ def _complex_hermitian_basis(n, field, normalize): Sij = _embed_complex_matrix(Eij) S.append(Sij) else: - # Beware, orthogonal but not normalized! The second one - # has a minus because it's conjugated. + # The second one has a minus because it's conjugated. Sij_real = _embed_complex_matrix(Eij + Eij.transpose()) S.append(Sij_real) Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose()) @@ -886,11 +920,7 @@ def _complex_hermitian_basis(n, field, normalize): # Since we embedded these, we can drop back to the "field" that we # started with instead of the complex extension "F". - S = [ s.change_ring(field) for s in S ] - if normalize: - S = [ s / _complex_hermitian_matrix_ip(s,s).sqrt() for s in S ] - - return tuple(S) + return tuple( s.change_ring(field) for s in S ) @@ -1206,10 +1236,6 @@ def _unembed_quaternion_matrix(M): return matrix(Q, n/4, elements) -# The usual inner product on R^n. -def _usual_ip(x,y): - return x.to_vector().inner_product(y.to_vector()) - # The inner product used for the real symmetric simple EJA. # We keep it as a separate function because e.g. the complex # algebra uses the same inner product, except divided by 2. @@ -1218,15 +1244,6 @@ def _matrix_ip(X,Y): Y_mat = Y.natural_representation() return (X_mat*Y_mat).trace() -def _real_symmetric_matrix_ip(X,Y): - return (X*Y).trace() - -def _complex_hermitian_matrix_ip(X,Y): - # This takes EMBEDDED matrices. - Xu = _unembed_complex_matrix(X) - Yu = _unembed_complex_matrix(Y) - # The trace need not be real; consider Xu = (i*I) and Yu = I. - return ((Xu*Yu).trace()).vector()[0] # real part, I guess class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra): """ @@ -1310,6 +1327,8 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra): """ def __init__(self, n, field=QQ, normalize_basis=True, **kwargs): + S = _real_symmetric_basis(n, field) + if n > 1 and normalize_basis: # We'll need sqrt(2) to normalize the basis, and this # winds up in the multiplication table, so the whole @@ -1319,8 +1338,12 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra): p = z**2 - 2 if p.is_irreducible(): field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt()) + S = [ s.change_ring(field) for s in S ] + self._basis_denormalizers = tuple( + self.__class__.natural_inner_product(s,s).sqrt() + for s in S ) + S = tuple( s/c for (s,c) in zip(S,self._basis_denormalizers) ) - S = _real_symmetric_basis(n, field, normalize_basis) Qs = _multiplication_table_from_matrix_basis(S) fdeja = super(RealSymmetricEJA, self) @@ -1330,10 +1353,6 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra): natural_basis=S, **kwargs) - def inner_product(self, x, y): - X = x.natural_representation() - Y = y.natural_representation() - return _real_symmetric_matrix_ip(X,Y) class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): @@ -1408,6 +1427,8 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): """ def __init__(self, n, field=QQ, normalize_basis=True, **kwargs): + S = _complex_hermitian_basis(n, field) + if n > 1 and normalize_basis: # We'll need sqrt(2) to normalize the basis, and this # winds up in the multiplication table, so the whole @@ -1417,8 +1438,12 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): p = z**2 - 2 if p.is_irreducible(): field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt()) + S = [ s.change_ring(field) for s in S ] + self._basis_denormalizers = tuple( + self.__class__.natural_inner_product(s,s).sqrt() + for s in S ) + S = tuple( s/c for (s,c) in zip(S,self._basis_denormalizers) ) - S = _complex_hermitian_basis(n, field, normalize_basis) Qs = _multiplication_table_from_matrix_basis(S) fdeja = super(ComplexHermitianEJA, self) @@ -1429,11 +1454,12 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): **kwargs) - def inner_product(self, x, y): - X = x.natural_representation() - Y = y.natural_representation() - return _complex_hermitian_matrix_ip(X,Y) - + @staticmethod + def natural_inner_product(X,Y): + Xu = _unembed_complex_matrix(X) + Yu = _unembed_complex_matrix(Y) + # The trace need not be real; consider Xu = (i*I) and Yu = I. + return ((Xu*Yu).trace()).vector()[0] # real part, I guess class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): """ @@ -1586,4 +1612,27 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra): return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs) def inner_product(self, x, y): - return _usual_ip(x,y) + """ + Faster to reimplement than to use natural representations. + + SETUP:: + + sage: from mjo.eja.eja_algebra import JordanSpinEJA + + TESTS: + + Ensure that this is the usual inner product for the algebras + over `R^n`:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: J = JordanSpinEJA(n) + sage: x = J.random_element() + sage: y = J.random_element() + sage: X = x.natural_representation() + sage: Y = y.natural_representation() + sage: x.inner_product(y) == J.__class__.natural_inner_product(X,Y) + True + + """ + return x.to_vector().inner_product(y.to_vector())