X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=106a0cddec06355a952e71d75699677b25dd7da9;hb=fc29add6cf1d9ff4e8a240b0f8f2ca6672d4ea57;hp=827d40214cac92cb6e1ee46eb17f4285a44398f6;hpb=a46720db62543983ab375654dee211ca844ac46c;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 827d402..106a0cd 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -172,6 +172,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): category = MagmaticAlgebras(field).FiniteDimensional() category = category.WithBasis().Unital().Commutative() + if n <= 1: + # All zero- and one-dimensional algebras are just the real + # numbers with (some positive multiples of) the usual + # multiplication as its Jordan and inner-product. + associative = True if associative is None: # We should figure it out. As with check_axioms, we have to do # this without the help of the _jordan_product_is_associative() @@ -1733,95 +1738,145 @@ class ConcreteEJA(FiniteDimensionalEJA): class MatrixEJA: @staticmethod - def jordan_product(X,Y): - return (X*Y + Y*X)/2 - - @staticmethod - def trace_inner_product(X,Y): - r""" - A trace inner-product for matrices that aren't embedded in the - reals. It takes MATRICES as arguments, not EJA elements. + def _denormalized_basis(A): """ - return (X*Y).trace().real() + Returns a basis for the space of complex Hermitian n-by-n matrices. -class RealEmbeddedMatrixEJA(MatrixEJA): - @staticmethod - def dimension_over_reals(): - r""" - The dimension of this matrix's base ring over the reals. + 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. - The reals are dimension one over themselves, obviously; that's - just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi` - have dimension two. Finally, the quaternions have dimension - four over the reals. + SETUP:: - This is used to determine the size of the matrix returned from - :meth:`real_embed`, among other things. - """ - raise NotImplementedError + sage: from mjo.hurwitz import (ComplexMatrixAlgebra, + ....: QuaternionMatrixAlgebra, + ....: OctonionMatrixAlgebra) + sage: from mjo.eja.eja_algebra import MatrixEJA - @classmethod - def real_embed(cls,M): - """ - Embed the matrix ``M`` into a space of real matrices. + TESTS:: - 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. + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: A = MatrixSpace(QQ, n) + sage: B = MatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B) + True - real_embed(M*N) = real_embed(M)*real_embed(N) + :: - """ - if M.ncols() != M.nrows(): - raise ValueError("the matrix 'M' must be square") - return M + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: A = ComplexMatrixAlgebra(n, scalars=QQ) + sage: B = MatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B) + True + :: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: A = QuaternionMatrixAlgebra(n, scalars=QQ) + sage: B = MatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B ) + True + + :: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: A = OctonionMatrixAlgebra(n, scalars=QQ) + sage: B = MatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B ) + True - @classmethod - def real_unembed(cls,M): - """ - The inverse of :meth:`real_embed`. """ - if M.ncols() != M.nrows(): - raise ValueError("the matrix 'M' must be square") - if not ZZ(M.nrows()).mod(cls.dimension_over_reals()).is_zero(): - raise ValueError("the matrix 'M' must be a real embedding") - return M + # These work for real MatrixSpace, whose monomials only have + # two coordinates (because the last one would always be "1"). + es = A.base_ring().gens() + gen = lambda A,m: A.monomial(m[:2]) + if hasattr(A, 'entry_algebra_gens'): + # We've got a MatrixAlgebra, and its monomials will have + # three coordinates. + es = A.entry_algebra_gens() + gen = lambda A,m: A.monomial(m) - @classmethod - def trace_inner_product(cls,X,Y): + basis = [] + for i in range(A.nrows()): + for j in range(i+1): + if i == j: + E_ii = gen(A, (i,j,es[0])) + basis.append(E_ii) + else: + for e in es: + E_ij = gen(A, (i,j,e)) + E_ij += E_ij.conjugate_transpose() + basis.append(E_ij) + + return tuple( basis ) + + @staticmethod + def jordan_product(X,Y): + return (X*Y + Y*X)/2 + + @staticmethod + def trace_inner_product(X,Y): r""" - Compute the trace inner-product of two real-embeddings. + A trace inner-product for matrices that aren't embedded in the + reals. It takes MATRICES as arguments, not EJA elements. SETUP:: - sage: from mjo.eja.eja_algebra import ComplexHermitianEJA + sage: from mjo.eja.eja_algebra import (RealSymmetricEJA, + ....: ComplexHermitianEJA, + ....: QuaternionHermitianEJA, + ....: OctonionHermitianEJA) EXAMPLES:: - sage: set_random_seed() - sage: J = ComplexHermitianEJA.random_instance() - sage: x,y = J.random_elements(2) - sage: Xe = x.to_matrix() - sage: Ye = y.to_matrix() - sage: X = J.real_unembed(Xe) - sage: Y = J.real_unembed(Ye) - sage: expected = (X*Y).trace().real() - sage: actual = J.trace_inner_product(Xe,Ye) - sage: actual == expected - True + sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 + + :: + + sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 + + :: + + sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 + + :: + + sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 """ - # This does in fact compute the real part of the trace. - # If we compute the trace of e.g. a complex matrix M, - # then we do so by adding up its diagonal entries -- - # call them z_1 through z_n. The real embedding of z_1 - # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace - # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth. - return (X*Y).trace()/cls.dimension_over_reals() + tr = (X*Y).trace() + if hasattr(tr, 'coefficient'): + # Works for octonions, and has to come first because they + # also have a "real()" method that doesn't return an + # element of the scalar ring. + return tr.coefficient(0) + elif hasattr(tr, 'coefficient_tuple'): + # Works for quaternions. + return tr.coefficient_tuple()[0] + + # Works for real and complex numbers. + return tr.real() + + class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): """ @@ -1888,38 +1943,6 @@ class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ - @classmethod - def _denormalized_basis(cls, n, field): - """ - Return a basis for the space of real symmetric n-by-n matrices. - - SETUP:: - - sage: from mjo.eja.eja_algebra import RealSymmetricEJA - - TESTS:: - - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ) - sage: all( M.is_symmetric() for M in B) - True - - """ - # The basis of symmetric matrices, as matrices, in their R^(n-by-n) - # coordinates. - S = [] - for i in range(n): - for j in range(i+1): - Eij = matrix(field, n, lambda k,l: k==i and l==j) - if i == j: - Sij = Eij - else: - Sij = Eij + Eij.transpose() - S.append(Sij) - return tuple(S) - - @staticmethod def _max_random_instance_size(): return 4 # Dimension 10 @@ -1937,175 +1960,22 @@ class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): # if the user passes check_axioms=True. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - associative = False - if n <= 1: - associative = True - - super().__init__(self._denormalized_basis(n,field), + A = MatrixSpace(field, n) + super().__init__(self._denormalized_basis(A), self.jordan_product, self.trace_inner_product, field=field, - associative=associative, **kwargs) # TODO: this could be factored out somehow, but is left here # because the MatrixEJA is not presently a subclass of the # FDEJA class that defines rank() and one(). self.rank.set_cache(n) - idV = self.matrix_space().one() - self.one.set_cache(self(idV)) - - - -class ComplexMatrixEJA(RealEmbeddedMatrixEJA): - # A manual dictionary-cache for the complex_extension() method, - # since apparently @classmethods can't also be @cached_methods. - _complex_extension = {} - - @classmethod - def complex_extension(cls,field): - r""" - The complex field that we embed/unembed, as an extension - of the given ``field``. - """ - if field in cls._complex_extension: - return cls._complex_extension[field] - - # Sage doesn't know how to adjoin the complex "i" (the root of - # x^2 + 1) to a field in a general way. Here, we just enumerate - # all of the cases that I have cared to support so far. - if field is AA: - # Sage doesn't know how to embed AA into QQbar, i.e. how - # to adjoin sqrt(-1) to AA. - F = QQbar - elif not field.is_exact(): - # RDF or RR - F = field.complex_field() - else: - # Works for QQ and... maybe some other fields. - R = PolynomialRing(field, 'z') - z = R.gen() - F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) - - cls._complex_extension[field] = F - return F - - @staticmethod - def dimension_over_reals(): - return 2 - - @classmethod - def real_embed(cls,M): - """ - Embed the n-by-n complex matrix ``M`` into the space of real - matrices of size 2n-by-2n via the map the sends each entry `z = a + - bi` to the block matrix ``[[a,b],[-b,a]]``. - - SETUP:: - - sage: from mjo.eja.eja_algebra import ComplexMatrixEJA + self.one.set_cache(self(A.one())) - EXAMPLES:: - - sage: F = QuadraticField(-1, 'I') - sage: x1 = F(4 - 2*i) - sage: x2 = F(1 + 2*i) - sage: x3 = F(-i) - sage: x4 = F(6) - sage: M = matrix(F,2,[[x1,x2],[x3,x4]]) - sage: ComplexMatrixEJA.real_embed(M) - [ 4 -2| 1 2] - [ 2 4|-2 1] - [-----+-----] - [ 0 -1| 6 0] - [ 1 0| 0 6] - TESTS: - Embedding is a homomorphism (isomorphism, in fact):: - - sage: set_random_seed() - sage: n = ZZ.random_element(3) - sage: F = QuadraticField(-1, 'I') - sage: X = random_matrix(F, n) - sage: Y = random_matrix(F, n) - sage: Xe = ComplexMatrixEJA.real_embed(X) - sage: Ye = ComplexMatrixEJA.real_embed(Y) - sage: XYe = ComplexMatrixEJA.real_embed(X*Y) - sage: Xe*Ye == XYe - True - - """ - super().real_embed(M) - n = M.nrows() - - # We don't need any adjoined elements... - field = M.base_ring().base_ring() - - blocks = [] - for z in M.list(): - a = z.real() - b = z.imag() - blocks.append(matrix(field, 2, [ [ a, b], - [-b, a] ])) - - return matrix.block(field, n, blocks) - - - @classmethod - def real_unembed(cls,M): - """ - The inverse of _embed_complex_matrix(). - - SETUP:: - - sage: from mjo.eja.eja_algebra import ComplexMatrixEJA - - EXAMPLES:: - - sage: A = matrix(QQ,[ [ 1, 2, 3, 4], - ....: [-2, 1, -4, 3], - ....: [ 9, 10, 11, 12], - ....: [-10, 9, -12, 11] ]) - sage: ComplexMatrixEJA.real_unembed(A) - [ 2*I + 1 4*I + 3] - [ 10*I + 9 12*I + 11] - - TESTS: - - Unembedding is the inverse of embedding:: - - sage: set_random_seed() - sage: F = QuadraticField(-1, 'I') - sage: M = random_matrix(F, 3) - sage: Me = ComplexMatrixEJA.real_embed(M) - sage: ComplexMatrixEJA.real_unembed(Me) == M - True - - """ - super().real_unembed(M) - n = ZZ(M.nrows()) - d = cls.dimension_over_reals() - F = cls.complex_extension(M.base_ring()) - i = F.gen() - - # Go top-left to bottom-right (reading order), converting every - # 2-by-2 block we see to a single complex element. - elements = [] - for k in range(n/d): - for j in range(n/d): - submat = M[d*k:d*k+d,d*j:d*j+d] - if submat[0,0] != submat[1,1]: - raise ValueError('bad on-diagonal submatrix') - if submat[0,1] != -submat[1,0]: - raise ValueError('bad off-diagonal submatrix') - z = submat[0,0] + submat[0,1]*i - elements.append(z) - - return matrix(F, n/d, elements) - - -class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, ComplexMatrixEJA): +class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, @@ -2162,90 +2032,24 @@ class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, ComplexMatrixEJA): Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ - - @classmethod - def _denormalized_basis(cls, 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 ComplexHermitianEJA - - TESTS:: - - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ) - sage: all( M.is_symmetric() for M in B) - True - - """ - R = PolynomialRing(ZZ, 'z') - z = R.gen() - F = ZZ.extension(z**2 + 1, 'I') - I = F.gen(1) - - # This is like the symmetric case, but we need to be careful: - # - # * We want conjugate-symmetry, not just symmetry. - # * The diagonal will (as a result) be real. - # - S = [] - Eij = matrix.zero(F,n) - for i in range(n): - for j in range(i+1): - # "build" E_ij - Eij[i,j] = 1 - if i == j: - Sij = cls.real_embed(Eij) - S.append(Sij) - else: - # The second one has a minus because it's conjugated. - Eij[j,i] = 1 # Eij = Eij + Eij.transpose() - Sij_real = cls.real_embed(Eij) - S.append(Sij_real) - # Eij = I*Eij - I*Eij.transpose() - Eij[i,j] = I - Eij[j,i] = -I - Sij_imag = cls.real_embed(Eij) - S.append(Sij_imag) - Eij[j,i] = 0 - # "erase" E_ij - Eij[i,j] = 0 - - # Since we embedded the entries, we can drop back to the - # desired real "field" instead of the extension "F". - return tuple( s.change_ring(field) for s in S ) - - def __init__(self, n, field=AA, **kwargs): # We know this is a valid EJA, but will double-check # if the user passes check_axioms=True. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - associative = False - if n <= 1: - associative = True - - super().__init__(self._denormalized_basis(n,field), + from mjo.hurwitz import ComplexMatrixAlgebra + A = ComplexMatrixAlgebra(n, scalars=field) + super().__init__(self._denormalized_basis(A), self.jordan_product, self.trace_inner_product, field=field, - associative=associative, **kwargs) + # TODO: this could be factored out somehow, but is left here # because the MatrixEJA is not presently a subclass of the # FDEJA class that defines rank() and one(). self.rank.set_cache(n) - idV = matrix.identity(ZZ, self.dimension_over_reals()*n) - self.one.set_cache(self(idV)) + self.one.set_cache(self(A.one())) @staticmethod def _max_random_instance_size(): @@ -2317,97 +2121,24 @@ class QuaternionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ - @classmethod - def _denormalized_basis(cls, 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 QuaternionHermitianEJA - - TESTS:: - - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ) - sage: all( M.is_hermitian() for M in B ) - True - - """ - from mjo.hurwitz import QuaternionMatrixAlgebra - A = QuaternionMatrixAlgebra(n, scalars=field) - es = A.entry_algebra_gens() - - basis = [] - for i in range(n): - for j in range(i+1): - if i == j: - E_ii = A.monomial( (i,j,es[0]) ) - basis.append(E_ii) - else: - for e in es: - E_ij = A.monomial( (i,j,e) ) - ec = e.conjugate() - # If the conjugate has a negative sign in front - # of it, (j,i,ec) won't be a monomial! - if (j,i,ec) in A.indices(): - E_ij += A.monomial( (j,i,ec) ) - else: - E_ij -= A.monomial( (j,i,-ec) ) - basis.append(E_ij) - - return tuple( basis ) - - - @staticmethod - def trace_inner_product(X,Y): - r""" - Overload the superclass method because the quaternions are weird - and we need to use ``coefficient_tuple()`` to get the realpart. - - SETUP:: - - sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA - - TESTS:: - - sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False) - sage: I = J.one().to_matrix() - sage: J.trace_inner_product(I, -I) - -2 - - """ - return (X*Y).trace().coefficient_tuple()[0] - def __init__(self, n, field=AA, **kwargs): # We know this is a valid EJA, but will double-check # if the user passes check_axioms=True. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - associative = False - if n <= 1: - associative = True - - super().__init__(self._denormalized_basis(n,field), + from mjo.hurwitz import QuaternionMatrixAlgebra + A = QuaternionMatrixAlgebra(n, scalars=field) + super().__init__(self._denormalized_basis(A), self.jordan_product, self.trace_inner_product, field=field, - associative=associative, **kwargs) # TODO: this could be factored out somehow, but is left here # because the MatrixEJA is not presently a subclass of the # FDEJA class that defines rank() and one(). self.rank.set_cache(n) - idV = self.matrix_space().one() - self.one.set_cache(self(idV)) + self.one.set_cache(self(A.one())) @staticmethod @@ -2534,7 +2265,9 @@ class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): # if the user passes check_axioms=True. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - super().__init__(self._denormalized_basis(n,field), + from mjo.hurwitz import OctonionMatrixAlgebra + A = OctonionMatrixAlgebra(n, scalars=field) + super().__init__(self._denormalized_basis(A), self.jordan_product, self.trace_inner_product, field=field, @@ -2544,72 +2277,7 @@ class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): # because the MatrixEJA is not presently a subclass of the # FDEJA class that defines rank() and one(). self.rank.set_cache(n) - idV = self.matrix_space().one() - self.one.set_cache(self(idV)) - - - @classmethod - def _denormalized_basis(cls, n, field): - """ - Returns a basis for the space of octonion Hermitian n-by-n - matrices. - - SETUP:: - - sage: from mjo.eja.eja_algebra import OctonionHermitianEJA - - EXAMPLES:: - - sage: B = OctonionHermitianEJA._denormalized_basis(3,QQ) - sage: all( M.is_hermitian() for M in B ) - True - sage: len(B) - 27 - - """ - from mjo.hurwitz import OctonionMatrixAlgebra - A = OctonionMatrixAlgebra(n, scalars=field) - es = A.entry_algebra_gens() - - basis = [] - for i in range(n): - for j in range(i+1): - if i == j: - E_ii = A.monomial( (i,j,es[0]) ) - basis.append(E_ii) - else: - for e in es: - E_ij = A.monomial( (i,j,e) ) - ec = e.conjugate() - # If the conjugate has a negative sign in front - # of it, (j,i,ec) won't be a monomial! - if (j,i,ec) in A.indices(): - E_ij += A.monomial( (j,i,ec) ) - else: - E_ij -= A.monomial( (j,i,-ec) ) - basis.append(E_ij) - - return tuple( basis ) - - @staticmethod - def trace_inner_product(X,Y): - r""" - The octonions don't know that the reals are embedded in them, - so we have to take the e0 component ourselves. - - SETUP:: - - sage: from mjo.eja.eja_algebra import OctonionHermitianEJA - - TESTS:: - - sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False) - sage: I = J.one().to_matrix() - sage: J.trace_inner_product(I, -I) - -2 - - """ - return (X*Y).trace().coefficient(0) + self.one.set_cache(self(A.one())) class AlbertEJA(OctonionHermitianEJA):