X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=79ccc7904830dc739e4af0a532c1bace2a801936;hb=9528af011cbb4d5e6a38ef972e0d14e7928d5eef;hp=a207250b4092f97e07d93a62c0ed7e23d9f9536d;hpb=7996413ef05ab2275c2cbb86494e30241904914b;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index a207250..79ccc79 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -5,6 +5,8 @@ are used in optimization, and have some additional nice methods beyond what can be supported in a general Jordan Algebra. """ +from itertools import izip, repeat + from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra from sage.categories.magmatic_algebras import MagmaticAlgebras from sage.combinat.free_module import CombinatorialFreeModule @@ -51,8 +53,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: set_random_seed() sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() + sage: x,y = J.random_elements(2) sage: x*y == y*x True @@ -60,9 +61,6 @@ 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() @@ -258,19 +256,6 @@ 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(): @@ -424,7 +409,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # assign a[r] goes out-of-bounds. a.append(1) # corresponds to x^r - return sum( a[k]*(t**k) for k in range(len(a)) ) + return sum( a[k]*(t**k) for k in xrange(len(a)) ) def inner_product(self, x, y): @@ -441,14 +426,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): EXAMPLES: - The inner product must satisfy its axiom for this algebra to truly - be a Euclidean Jordan Algebra:: + Our inner product satisfies the Jordan axiom, which is also + referred to as "associativity" for a symmetric bilinear form:: sage: set_random_seed() sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() + sage: x,y,z = J.random_elements(3) sage: (x*y).inner_product(z) == y.inner_product(x*z) True @@ -508,7 +491,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): """ M = list(self._multiplication_table) # copy - for i in range(len(M)): + for i in xrange(len(M)): # M had better be "square" M[i] = [self.monomial(i)] + M[i] M = [["*"] + list(self.gens())] + M @@ -657,6 +640,25 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): s = super(FiniteDimensionalEuclideanJordanAlgebra, self) return s.random_element() + def random_elements(self, count): + """ + Return ``count`` random elements as a tuple. + + SETUP:: + + sage: from mjo.eja.eja_algebra import JordanSpinEJA + + EXAMPLES:: + + sage: J = JordanSpinEJA(3) + sage: x,y,z = J.random_elements(3) + sage: all( [ x in J, y in J, z in J ]) + True + sage: len( J.random_elements(10) ) == 10 + True + + """ + return tuple( self.random_element() for idx in xrange(count) ) @classmethod def random_instance(cls, field=QQ, **kwargs): @@ -676,7 +678,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # not worry about it. raise NotImplementedError - n = ZZ.random_element(1, cls._max_test_case_size()) + n = ZZ.random_element(cls._max_test_case_size()) + 1 return cls(n, field, **kwargs) @@ -713,18 +715,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): The rank of the `n`-by-`n` Hermitian real, complex, or quaternion matrices is `n`:: - sage: RealSymmetricEJA(2).rank() - 2 - sage: ComplexHermitianEJA(2).rank() - 2 + sage: RealSymmetricEJA(4).rank() + 4 + sage: ComplexHermitianEJA(3).rank() + 3 sage: QuaternionHermitianEJA(2).rank() 2 - sage: RealSymmetricEJA(5).rank() - 5 - sage: ComplexHermitianEJA(5).rank() - 5 - sage: QuaternionHermitianEJA(5).rank() - 5 TESTS: @@ -800,21 +796,11 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: RealCartesianProductEJA(3, prefix='r').gens() (r0, r1, r2) - Our inner product satisfies the Jordan axiom:: - - sage: set_random_seed() - sage: J = RealCartesianProductEJA.random_instance() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - sage: (x*y).inner_product(z) == y.inner_product(x*z) - True - """ def __init__(self, n, field=QQ, **kwargs): V = VectorSpace(field, n) - mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ] - for i in range(n) ] + mult_table = [ [ V.gen(i)*(i == j) for j in xrange(n) ] + for i in xrange(n) ] fdeja = super(RealCartesianProductEJA, self) return fdeja.__init__(field, mult_table, rank=n, **kwargs) @@ -834,8 +820,7 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: set_random_seed() sage: J = RealCartesianProductEJA.random_instance() - sage: x = J.random_element() - sage: y = J.random_element() + sage: x,y = J.random_elements(2) sage: X = x.natural_representation() sage: Y = y.natural_representation() sage: x.inner_product(y) == J.natural_inner_product(X,Y) @@ -898,7 +883,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): def _max_test_case_size(): # Play it safe, since this will be squared and the underlying # field can have dimension 4 (quaternions) too. - return 3 + return 2 @classmethod def _denormalized_basis(cls, n, field): @@ -907,6 +892,9 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): def __init__(self, n, field=QQ, normalize_basis=True, **kwargs): S = self._denormalized_basis(n, field) + # Used in this class's fast _charpoly_coeff() override. + self._basis_normalizers = None + 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 @@ -931,6 +919,30 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): **kwargs) + @cached_method + def _charpoly_coeff(self, i): + """ + Override the parent method with something that tries to compute + over a faster (non-extension) field. + """ + if self._basis_normalizers is None: + # We didn't normalize, so assume that the basis we started + # with had entries in a nice field. + return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i) + else: + # If we didn't unembed first, this number would be wrong + # by a power-of-two factor for complex/quaternion matrices. + n = self.real_unembed(self.natural_basis_space().zero()).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 = izip(x.base_ring().gens(), self._basis_normalizers) + substitutions = { v: v*c for (v,c) in pairs } + return p.subs(substitutions) + + @staticmethod def multiplication_table_from_matrix_basis(basis): """ @@ -952,9 +964,9 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): V = VectorSpace(field, dimension**2) W = V.span_of_basis( _mat2vec(s) for s in basis ) n = len(basis) - mult_table = [[W.zero() for j in range(n)] for i in range(n)] - for i in range(n): - for j in range(n): + mult_table = [[W.zero() for j in xrange(n)] for i in xrange(n)] + for i in xrange(n): + for j in xrange(n): mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2 mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry)) @@ -1006,7 +1018,33 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): return tr.coefficient_tuple()[0] -class RealSymmetricEJA(MatrixEuclideanJordanAlgebra): +class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): + @staticmethod + def real_embed(M): + """ + Embed the matrix ``M`` into a space of real matrices. + + The matrix ``M`` can have entries in any field at the moment: + the real numbers, complex numbers, or quaternions. And although + they are not a field, we can probably support octonions at some + point, too. This function returns a real matrix that "acts like" + the original with respect to matrix multiplication; i.e. + + real_embed(M*N) = real_embed(M)*real_embed(N) + + """ + return M + + + @staticmethod + def real_unembed(M): + """ + The inverse of :meth:`real_embed`. + """ + return M + + +class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -1042,8 +1080,7 @@ class RealSymmetricEJA(MatrixEuclideanJordanAlgebra): sage: set_random_seed() sage: J = RealSymmetricEJA.random_instance() - sage: x = J.random_element() - sage: y = J.random_element() + sage: x,y = J.random_elements(2) sage: actual = (x*y).natural_representation() sage: X = x.natural_representation() sage: Y = y.natural_representation() @@ -1058,16 +1095,6 @@ class RealSymmetricEJA(MatrixEuclideanJordanAlgebra): sage: RealSymmetricEJA(3, prefix='q').gens() (q0, q1, q2, q3, q4, q5) - Our inner product satisfies the Jordan axiom:: - - sage: set_random_seed() - sage: J = RealSymmetricEJA.random_instance() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - sage: (x*y).inner_product(z) == y.inner_product(x*z) - True - Our natural basis is normalized with respect to the natural inner product unless we specify otherwise:: @@ -1122,31 +1149,7 @@ class RealSymmetricEJA(MatrixEuclideanJordanAlgebra): @staticmethod def _max_test_case_size(): - return 5 # Dimension 10 - - @staticmethod - def real_embed(M): - """ - Embed the matrix ``M`` into a space of real matrices. - - The matrix ``M`` can have entries in any field at the moment: - the real numbers, complex numbers, or quaternions. And although - they are not a field, we can probably support octonions at some - point, too. This function returns a real matrix that "acts like" - the original with respect to matrix multiplication; i.e. - - real_embed(M*N) = real_embed(M)*real_embed(N) - - """ - return M - - - @staticmethod - def real_unembed(M): - """ - The inverse of :meth:`real_embed`. - """ - return M + return 4 # Dimension 10 @@ -1269,6 +1272,37 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): return matrix(F, n/2, elements) + @classmethod + def natural_inner_product(cls,X,Y): + """ + Compute a natural inner product in this algebra directly from + its real embedding. + + SETUP:: + + sage: from mjo.eja.eja_algebra import ComplexHermitianEJA + + TESTS: + + This gives the same answer as the slow, default method implemented + in :class:`MatrixEuclideanJordanAlgebra`:: + + sage: set_random_seed() + sage: J = ComplexHermitianEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: Xe = x.natural_representation() + sage: Ye = y.natural_representation() + sage: X = ComplexHermitianEJA.real_unembed(Xe) + sage: Y = ComplexHermitianEJA.real_unembed(Ye) + sage: expected = (X*Y).trace().vector()[0] + sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye) + sage: actual == expected + True + + """ + return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2 + + class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of complex Hermitian n-by-n @@ -1295,8 +1329,7 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() - sage: x = J.random_element() - sage: y = J.random_element() + sage: x,y = J.random_elements(2) sage: actual = (x*y).natural_representation() sage: X = x.natural_representation() sage: Y = y.natural_representation() @@ -1311,16 +1344,6 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): sage: ComplexHermitianEJA(2, prefix='z').gens() (z0, z1, z2, z3) - Our inner product satisfies the Jordan axiom:: - - sage: set_random_seed() - sage: J = ComplexHermitianEJA.random_instance() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - sage: (x*y).inner_product(z) == y.inner_product(x*z) - True - Our natural basis is normalized with respect to the natural inner product unless we specify otherwise:: @@ -1528,6 +1551,36 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): return matrix(Q, n/4, elements) + @classmethod + def natural_inner_product(cls,X,Y): + """ + Compute a natural inner product in this algebra directly from + its real embedding. + + SETUP:: + + sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA + + TESTS: + + This gives the same answer as the slow, default method implemented + in :class:`MatrixEuclideanJordanAlgebra`:: + + sage: set_random_seed() + sage: J = QuaternionHermitianEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: Xe = x.natural_representation() + sage: Ye = y.natural_representation() + sage: X = QuaternionHermitianEJA.real_unembed(Xe) + sage: Y = QuaternionHermitianEJA.real_unembed(Ye) + sage: expected = (X*Y).trace().coefficient_tuple()[0] + sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye) + sage: actual == expected + True + + """ + return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4 + class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): """ @@ -1555,8 +1608,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): sage: set_random_seed() sage: J = QuaternionHermitianEJA.random_instance() - sage: x = J.random_element() - sage: y = J.random_element() + sage: x,y = J.random_elements(2) sage: actual = (x*y).natural_representation() sage: X = x.natural_representation() sage: Y = y.natural_representation() @@ -1571,16 +1623,6 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): sage: QuaternionHermitianEJA(2, prefix='a').gens() (a0, a1, a2, a3, a4, a5) - Our inner product satisfies the Jordan axiom:: - - sage: set_random_seed() - sage: J = QuaternionHermitianEJA.random_instance() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - sage: (x*y).inner_product(z) == y.inner_product(x*z) - True - Our natural basis is normalized with respect to the natural inner product unless we specify otherwise:: @@ -1692,22 +1734,12 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: JordanSpinEJA(2, prefix='B').gens() (B0, B1) - Our inner product satisfies the Jordan axiom:: - - sage: set_random_seed() - sage: J = JordanSpinEJA.random_instance() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - sage: (x*y).inner_product(z) == y.inner_product(x*z) - True - """ def __init__(self, n, field=QQ, **kwargs): V = VectorSpace(field, n) - mult_table = [[V.zero() for j in range(n)] for i in range(n)] - for i in range(n): - for j in range(n): + mult_table = [[V.zero() for j in xrange(n)] for i in xrange(n)] + for i in xrange(n): + for j in xrange(n): x = V.gen(i) y = V.gen(j) x0 = x[0] @@ -1741,8 +1773,7 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: set_random_seed() sage: J = JordanSpinEJA.random_instance() - sage: x = J.random_element() - sage: y = J.random_element() + sage: x,y = J.random_elements(2) sage: X = x.natural_representation() sage: Y = y.natural_representation() sage: x.inner_product(y) == J.natural_inner_product(X,Y)