From: Michael Orlitzky Date: Tue, 9 Mar 2021 04:58:39 +0000 (-0500) Subject: eja: refactor all matrix classes upwards (note: everything broken). X-Git-Url: http://gitweb.michael.orlitzky.com/?p=sage.d.git;a=commitdiff_plain;h=fc29add6cf1d9ff4e8a240b0f8f2ca6672d4ea57 eja: refactor all matrix classes upwards (note: everything broken). --- diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 3027075..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() @@ -1732,6 +1737,86 @@ class ConcreteEJA(FiniteDimensionalEJA): class MatrixEJA: + @staticmethod + def _denormalized_basis(A): + """ + 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.hurwitz import (ComplexMatrixAlgebra, + ....: QuaternionMatrixAlgebra, + ....: OctonionMatrixAlgebra) + sage: from mjo.eja.eja_algebra import MatrixEJA + + TESTS:: + + 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 + + :: + + 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 + + """ + # 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) + + 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 @@ -1741,8 +1826,56 @@ class MatrixEJA: r""" 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 (RealSymmetricEJA, + ....: ComplexHermitianEJA, + ....: QuaternionHermitianEJA, + ....: OctonionHermitianEJA) + + EXAMPLES:: + + 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 + """ - return (X*Y).trace().real() + 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): @@ -1810,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 @@ -1859,23 +1960,18 @@ 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)) + self.one.set_cache(self(A.one())) @@ -1936,93 +2032,24 @@ class ComplexHermitianEJA(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 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,QQ) - sage: all( M.is_hermitian() for M in B) - True - - """ - from mjo.hurwitz import ComplexMatrixAlgebra - A = ComplexMatrixAlgebra(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""" - SETUP:: - - sage: from mjo.eja.eja_algebra import ComplexHermitianEJA - - TESTS:: - - sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) - sage: I = J.one().to_matrix() - sage: J.trace_inner_product(I, -I) - -2 - - """ - return (X*Y).trace().real() - 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 = self.matrix_space().one() - self.one.set_cache(self(idV)) + self.one.set_cache(self(A.one())) @staticmethod def _max_random_instance_size(): @@ -2094,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 @@ -2311,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, @@ -2321,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): diff --git a/mjo/hurwitz.py b/mjo/hurwitz.py index ccc8219..1f7c9dc 100644 --- a/mjo/hurwitz.py +++ b/mjo/hurwitz.py @@ -306,6 +306,31 @@ class Octonions(CombinatorialFreeModule): class HurwitzMatrixAlgebraElement(MatrixAlgebraElement): + def conjugate_transpose(self): + r""" + Return the conjugate-transpose of this matrix. + + SETUP:: + + sage: from mjo.hurwitz import HurwitzMatrixAlgebra + + EXAMPLES:: + + sage: A = HurwitzMatrixAlgebra(2, QQbar, ZZ) + sage: M = A([ [ I, 2*I], + ....: [ 3*I, 4*I] ]) + +------+------+ + | -1*I | -3*I | + +------+------+ + | -2*I | -4*I | + +------+------+ + + """ + entries = [ [ self[j,i].conjugate() + for j in range(self.ncols())] + for i in range(self.nrows()) ] + return self.parent()._element_constructor_(entries) + def is_hermitian(self): r""" @@ -322,6 +347,8 @@ class HurwitzMatrixAlgebraElement(MatrixAlgebraElement): True """ + # A tiny bit faster than checking equality with the conjugate + # transpose. return all( self[i,j] == self[j,i].conjugate() for i in range(self.nrows()) for j in range(self.ncols()) )