X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=4947fe82bb981401a09a039306516b4b0631113a;hb=f911ca107b6ad2a92547a255311ef7f16978feac;hp=07c49473a044411374583b60a15127096c355c2c;hpb=4e1d889fa70b660e789f31adb13f449c9f3a2fb1;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 07c4947..4947fe8 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -15,9 +15,10 @@ from sage.misc.prandom import choice from sage.misc.table import table from sage.modules.free_module import FreeModule, VectorSpace from sage.rings.integer_ring import ZZ -from sage.rings.number_field.number_field import QuadraticField +from sage.rings.number_field.number_field import NumberField, QuadraticField from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.rational_field import QQ +from sage.rings.real_lazy import CLF, RLF from sage.structure.element import is_Matrix from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement @@ -59,6 +60,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): self._rank = rank self._natural_basis = natural_basis + # TODO: HACK for the charpoly.. needs redesign badly. + self._basis_normalizers = None + if category is None: category = MagmaticAlgebras(field).FiniteDimensional() category = category.WithBasis().Unital() @@ -135,13 +139,17 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return self.zero() natural_basis = self.natural_basis() - if elt not in natural_basis[0].matrix_space(): + basis_space = natural_basis[0].matrix_space() + if elt not in basis_space: raise ValueError("not a naturally-represented algebra element") - # Thanks for nothing! Matrix spaces aren't vector - # spaces in Sage, so we have to figure out its - # natural-basis coordinates ourselves. - V = VectorSpace(elt.base_ring(), elt.nrows()*elt.ncols()) + # Thanks for nothing! Matrix spaces aren't vector spaces in + # Sage, so we have to figure out its natural-basis coordinates + # ourselves. We use the basis space's ring instead of the + # element's ring because the basis space might be an algebraic + # closure whereas the base ring of the 3-by-3 identity matrix + # could be QQ instead of QQbar. + V = VectorSpace(basis_space.base_ring(), elt.nrows()*elt.ncols()) W = V.span_of_basis( _mat2vec(s) for s in natural_basis ) coords = W.coordinate_vector(_mat2vec(elt)) return self.from_vector(coords) @@ -219,6 +227,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return V.span_of_basis(b) + @cached_method def _charpoly_coeff(self, i): """ @@ -229,6 +238,19 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): store the trace/determinant (a_{r-1} and a_{0} respectively) separate from the entire characteristic polynomial. """ + if self._basis_normalizers is not None: + # Must be a matrix class? + # WARNING/TODO: this whole mess is mis-designed. + n = self.natural_basis_space().nrows() + field = self.base_ring().base_ring() # yeeeeaaaahhh + J = self.__class__(n, field, False) + (_,x,_,_) = J._charpoly_matrix_system() + p = J._charpoly_coeff(i) + # p might be missing some vars, have to substitute "optionally" + pairs = zip(x.base_ring().gens(), self._basis_normalizers) + substitutions = { v: v*c for (v,c) in pairs } + return p.subs(substitutions) + (A_of_x, x, xr, detA) = self._charpoly_matrix_system() R = A_of_x.base_ring() if i >= self.rank(): @@ -411,9 +433,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): @@ -498,8 +520,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Finite family {0: e0, 1: e1, 2: e2} sage: J.natural_basis() ( - [1 0] [0 1] [0 0] - [0 0], [1 0], [0 1] + [1 0] [ 0 1/2*sqrt2] [0 0] + [0 0], [1/2*sqrt2 0], [0 1] ) :: @@ -532,6 +554,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): """ @@ -674,7 +710,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = RealSymmetricEJA(2) sage: J.vector_space() - Vector space of dimension 3 over Rational Field + Vector space of dimension 3 over... """ return self.zero().to_vector().parent().ambient_vector_space() @@ -743,7 +779,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(): @@ -800,6 +859,19 @@ def random_eja(): def _real_symmetric_basis(n, field): """ Return a basis for the space of real symmetric n-by-n matrices. + + SETUP:: + + sage: from mjo.eja.eja_algebra import _real_symmetric_basis + + TESTS:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: B = _real_symmetric_basis(n, QQ) + sage: all( M.is_symmetric() for M in B) + True + """ # The basis of symmetric matrices, as matrices, in their R^(n-by-n) # coordinates. @@ -810,7 +882,6 @@ def _real_symmetric_basis(n, field): if i == j: Sij = Eij else: - # Beware, orthogonal but not normalized! Sij = Eij + Eij.transpose() S.append(Sij) return tuple(S) @@ -820,6 +891,12 @@ def _complex_hermitian_basis(n, field): """ Returns a basis for the space of complex Hermitian n-by-n matrices. + Why do we embed these? Basically, because all of numerical linear + algebra assumes that you're working with vectors consisting of `n` + entries from a field and scalars from the same field. There's no way + to tell SageMath that (for example) the vectors contain complex + numbers, while the scalar field is real. + SETUP:: sage: from mjo.eja.eja_algebra import _complex_hermitian_basis @@ -828,11 +905,15 @@ def _complex_hermitian_basis(n, field): sage: set_random_seed() sage: n = ZZ.random_element(1,5) - sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) ) + sage: field = QuadraticField(2, 'sqrt2') + sage: B = _complex_hermitian_basis(n, field) + sage: all( M.is_symmetric() for M in B) True """ - F = QuadraticField(-1, 'I') + R = PolynomialRing(field, 'z') + z = R.gen() + F = NumberField(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) I = F.gen() # This is like the symmetric case, but we need to be careful: @@ -843,24 +924,33 @@ def _complex_hermitian_basis(n, field): S = [] for i in xrange(n): for j in xrange(i+1): - Eij = matrix(field, n, lambda k,l: k==i and l==j) + Eij = matrix(F, n, lambda k,l: k==i and l==j) if i == j: 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()) S.append(Sij_imag) - return tuple(S) + + # Since we embedded these, we can drop back to the "field" that we + # started with instead of the complex extension "F". + return tuple( s.change_ring(field) for s in S ) + def _quaternion_hermitian_basis(n, field): """ Returns a basis for the space of quaternion Hermitian n-by-n matrices. + Why do we embed these? Basically, because all of numerical linear + algebra assumes that you're working with vectors consisting of `n` + entries from a field and scalars from the same field. There's no way + to tell SageMath that (for example) the vectors contain complex + numbers, while the scalar field is real. + SETUP:: sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis @@ -869,7 +959,8 @@ def _quaternion_hermitian_basis(n, field): sage: set_random_seed() sage: n = ZZ.random_element(1,5) - sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) ) + sage: B = _quaternion_hermitian_basis(n, QQ) + sage: all( M.is_symmetric() for M in B ) True """ @@ -945,7 +1036,7 @@ def _embed_complex_matrix(M): EXAMPLES:: - sage: F = QuadraticField(-1,'i') + sage: F = QuadraticField(-1, 'i') sage: x1 = F(4 - 2*i) sage: x2 = F(1 + 2*i) sage: x3 = F(-i) @@ -979,8 +1070,8 @@ def _embed_complex_matrix(M): field = M.base_ring() blocks = [] for z in M.list(): - a = z.real() - b = z.imag() + a = z.vector()[0] # real part, I guess + b = z.vector()[1] # imag part, I guess blocks.append(matrix(field, 2, [[a,b],[-b,a]])) # We can drop the imaginaries here. @@ -1023,7 +1114,10 @@ def _unembed_complex_matrix(M): if not n.mod(2).is_zero(): raise ValueError("the matrix 'M' must be a complex embedding") - F = QuadraticField(-1, 'i') + field = M.base_ring() # This should already have sqrt2 + R = PolynomialRing(field, 'z') + z = R.gen() + F = NumberField(z**2 + 1,'i', embedding=CLF(-1).sqrt()) i = F.gen() # Go top-left to bottom-right (reading order), converting every @@ -1159,10 +1253,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. @@ -1189,7 +1279,7 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: e0*e0 e0 sage: e1*e1 - e0 + e2 + 1/2*e0 + 1/2*e2 sage: e2*e2 e2 @@ -1235,9 +1325,42 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: (x*y).inner_product(z) == y.inner_product(x*z) True + Our basis is normalized with respect to the natural inner product:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: J = RealSymmetricEJA(n) + sage: all( b.norm() == 1 for b in J.gens() ) + True + + Left-multiplication operators are symmetric because they satisfy + the Jordan axiom:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: x = RealSymmetricEJA(n).random_element() + sage: x.operator().matrix().is_symmetric() + True + """ - def __init__(self, n, field=QQ, **kwargs): + 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 + # algebra needs to be over the field extension. + R = PolynomialRing(field, 'z') + z = R.gen() + 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_normalizers = 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_normalizers) ) + Qs = _multiplication_table_from_matrix_basis(S) fdeja = super(RealSymmetricEJA, self) @@ -1247,8 +1370,6 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra): natural_basis=S, **kwargs) - def inner_product(self, x, y): - return _matrix_ip(x,y) class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): @@ -1304,9 +1425,42 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: (x*y).inner_product(z) == y.inner_product(x*z) True + Our basis is normalized with respect to the natural inner product:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,4) + sage: J = ComplexHermitianEJA(n) + sage: all( b.norm() == 1 for b in J.gens() ) + True + + Left-multiplication operators are symmetric because they satisfy + the Jordan axiom:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: x = ComplexHermitianEJA(n).random_element() + sage: x.operator().matrix().is_symmetric() + True + """ - def __init__(self, n, field=QQ, **kwargs): + 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 + # algebra needs to be over the field extension. + R = PolynomialRing(field, 'z') + z = R.gen() + 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_normalizers = 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_normalizers) ) + Qs = _multiplication_table_from_matrix_basis(S) fdeja = super(ComplexHermitianEJA, self) @@ -1317,16 +1471,12 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): **kwargs) - def inner_product(self, x, y): - # Since a+bi on the diagonal is represented as - # - # a + bi = [ a b ] - # [ -b a ], - # - # we'll double-count the "a" entries if we take the trace of - # the embedding. - return _matrix_ip(x,y)/2 - + @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): """ @@ -1382,8 +1532,8 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): True """ - def __init__(self, n, field=QQ, **kwargs): - S = _quaternion_hermitian_basis(n, field) + def __init__(self, n, field=QQ, normalize_basis=True, **kwargs): + S = _quaternion_hermitian_basis(n, field, normalize_basis) Qs = _multiplication_table_from_matrix_basis(S) fdeja = super(QuaternionHermitianEJA, self) @@ -1479,4 +1629,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())