From a46720db62543983ab375654dee211ca844ac46c Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 8 Mar 2021 16:00:31 -0500 Subject: [PATCH] eja: don't real-embed quaternion matrices. --- mjo/all.py | 2 +- mjo/eja/eja_algebra.py | 271 ++++++++--------------------------------- mjo/eja/eja_utils.py | 18 ++- 3 files changed, 71 insertions(+), 220 deletions(-) diff --git a/mjo/all.py b/mjo/all.py index f802767..78550cc 100644 --- a/mjo/all.py +++ b/mjo/all.py @@ -8,7 +8,7 @@ from mjo.cone.all import * from mjo.eja.all import * from mjo.interpolation import * from mjo.misc import * -from mjo.octonions import Octonions +from mjo.hurwitz import Octonions from mjo.orthogonal_polynomials import * from mjo.polynomial import * from mjo.symbol_sequence import * diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 2ecf8cd..827d402 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -1019,7 +1019,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Full MatrixSpace of 4 by 4 dense matrices over Rational Field sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False) sage: J.matrix_space() - Full MatrixSpace of 4 by 4 dense matrices over Rational Field + Module of 1 by 1 matrices with entries in Quaternion + Algebra (-1, -1) with base ring Rational Field over + the scalar ring Rational Field """ if self.is_trivial(): @@ -1796,8 +1798,7 @@ class RealEmbeddedMatrixEJA(MatrixEJA): SETUP:: - sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, - ....: QuaternionHermitianEJA) + sage: from mjo.eja.eja_algebra import ComplexHermitianEJA EXAMPLES:: @@ -1813,20 +1814,6 @@ class RealEmbeddedMatrixEJA(MatrixEJA): sage: actual == expected True - :: - - sage: set_random_seed() - sage: J = QuaternionHermitianEJA.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().coefficient_tuple()[0] - sage: actual = J.trace_inner_product(Xe,Ye) - sage: actual == expected - True - """ # This does in fact compute the real part of the trace. # If we compute the trace of e.g. a complex matrix M, @@ -2272,157 +2259,8 @@ class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, ComplexMatrixEJA): n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, **kwargs) -class QuaternionMatrixEJA(RealEmbeddedMatrixEJA): - - # A manual dictionary-cache for the quaternion_extension() method, - # since apparently @classmethods can't also be @cached_methods. - _quaternion_extension = {} - - @classmethod - def quaternion_extension(cls,field): - r""" - The quaternion field that we embed/unembed, as an extension - of the given ``field``. - """ - if field in cls._quaternion_extension: - return cls._quaternion_extension[field] - - Q = QuaternionAlgebra(field,-1,-1) - - cls._quaternion_extension[field] = Q - return Q - - @staticmethod - def dimension_over_reals(): - return 4 - - @classmethod - def real_embed(cls,M): - """ - Embed the n-by-n quaternion matrix ``M`` into the space of real - matrices of size 4n-by-4n by first sending each quaternion entry `z - = a + bi + cj + dk` to the block-complex matrix ``[[a + bi, - c+di],[-c + di, a-bi]]`, and then embedding those into a real - matrix. - - SETUP:: - - sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA - - EXAMPLES:: - - sage: Q = QuaternionAlgebra(QQ,-1,-1) - sage: i,j,k = Q.gens() - sage: x = 1 + 2*i + 3*j + 4*k - sage: M = matrix(Q, 1, [[x]]) - sage: QuaternionMatrixEJA.real_embed(M) - [ 1 2 3 4] - [-2 1 -4 3] - [-3 4 1 -2] - [-4 -3 2 1] - - Embedding is a homomorphism (isomorphism, in fact):: - - sage: set_random_seed() - sage: n = ZZ.random_element(2) - sage: Q = QuaternionAlgebra(QQ,-1,-1) - sage: X = random_matrix(Q, n) - sage: Y = random_matrix(Q, n) - sage: Xe = QuaternionMatrixEJA.real_embed(X) - sage: Ye = QuaternionMatrixEJA.real_embed(Y) - sage: XYe = QuaternionMatrixEJA.real_embed(X*Y) - sage: Xe*Ye == XYe - True - - """ - super().real_embed(M) - quaternions = M.base_ring() - n = M.nrows() - - F = QuadraticField(-1, 'I') - i = F.gen() - - blocks = [] - for z in M.list(): - t = z.coefficient_tuple() - a = t[0] - b = t[1] - c = t[2] - d = t[3] - cplxM = matrix(F, 2, [[ a + b*i, c + d*i], - [-c + d*i, a - b*i]]) - realM = ComplexMatrixEJA.real_embed(cplxM) - blocks.append(realM) - - # We should have real entries by now, so use the realest field - # we've got for the return value. - return matrix.block(quaternions.base_ring(), n, blocks) - - - - @classmethod - def real_unembed(cls,M): - """ - The inverse of _embed_quaternion_matrix(). - - SETUP:: - - sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA - - EXAMPLES:: - - sage: M = matrix(QQ, [[ 1, 2, 3, 4], - ....: [-2, 1, -4, 3], - ....: [-3, 4, 1, -2], - ....: [-4, -3, 2, 1]]) - sage: QuaternionMatrixEJA.real_unembed(M) - [1 + 2*i + 3*j + 4*k] - - TESTS: - - Unembedding is the inverse of embedding:: - - sage: set_random_seed() - sage: Q = QuaternionAlgebra(QQ, -1, -1) - sage: M = random_matrix(Q, 3) - sage: Me = QuaternionMatrixEJA.real_embed(M) - sage: QuaternionMatrixEJA.real_unembed(Me) == M - True - - """ - super().real_unembed(M) - n = ZZ(M.nrows()) - d = cls.dimension_over_reals() - - # Use the base ring of the matrix to ensure that its entries can be - # multiplied by elements of the quaternion algebra. - Q = cls.quaternion_extension(M.base_ring()) - i,j,k = Q.gens() - - # Go top-left to bottom-right (reading order), converting every - # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1 - # quaternion block. - elements = [] - for l in range(n/d): - for m in range(n/d): - submat = ComplexMatrixEJA.real_unembed( - M[d*l:d*l+d,d*m:d*m+d] ) - if submat[0,0] != submat[1,1].conjugate(): - raise ValueError('bad on-diagonal submatrix') - if submat[0,1] != -submat[1,0].conjugate(): - raise ValueError('bad off-diagonal submatrix') - z = submat[0,0].real() - z += submat[0,0].imag()*i - z += submat[0,1].real()*j - z += submat[0,1].imag()*k - elements.append(z) - - return matrix(Q, n/d, elements) - -class QuaternionHermitianEJA(RationalBasisEJA, - ConcreteEJA, - QuaternionMatrixEJA): +class QuaternionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): r""" The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the @@ -2498,59 +2336,56 @@ class QuaternionHermitianEJA(RationalBasisEJA, sage: set_random_seed() sage: n = ZZ.random_element(1,5) - sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ) - sage: all( M.is_symmetric() for M in B ) + sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ) + sage: all( M.is_hermitian() for M in B ) True """ - Q = QuaternionAlgebra(QQ,-1,-1) - I,J,K = Q.gens() + from mjo.hurwitz import QuaternionMatrixAlgebra + A = QuaternionMatrixAlgebra(n, scalars=field) + es = A.entry_algebra_gens() - # 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(Q,n) + basis = [] 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) + E_ii = A.monomial( (i,j,es[0]) ) + basis.append(E_ii) else: - # The second, third, and fourth ones have a minus - # because they're conjugated. - # Eij = Eij + Eij.transpose() - Eij[j,i] = 1 - Sij_real = cls.real_embed(Eij) - S.append(Sij_real) - # Eij = I*(Eij - Eij.transpose()) - Eij[i,j] = I - Eij[j,i] = -I - Sij_I = cls.real_embed(Eij) - S.append(Sij_I) - # Eij = J*(Eij - Eij.transpose()) - Eij[i,j] = J - Eij[j,i] = -J - Sij_J = cls.real_embed(Eij) - S.append(Sij_J) - # Eij = K*(Eij - Eij.transpose()) - Eij[i,j] = K - Eij[j,i] = -K - Sij_K = cls.real_embed(Eij) - S.append(Sij_K) - Eij[j,i] = 0 - # "erase" E_ij - Eij[i,j] = 0 + 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) - # Since we embedded the entries, we can drop back to the - # desired real "field" instead of the quaternion algebra "Q". - return tuple( s.change_ring(field) for s in S ) + 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. @@ -2571,7 +2406,7 @@ class QuaternionHermitianEJA(RationalBasisEJA, # 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) + idV = self.matrix_space().one() self.one.set_cache(self(idV)) @@ -2733,25 +2568,25 @@ class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): """ from mjo.hurwitz import OctonionMatrixAlgebra - MS = OctonionMatrixAlgebra(n, scalars=field) - es = MS.entry_algebra().gens() + 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 = MS.monomial( (i,j,es[0]) ) + E_ii = A.monomial( (i,j,es[0]) ) basis.append(E_ii) else: for e in es: - E_ij = MS.monomial( (i,j,e) ) + 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 MS.indices(): - E_ij += MS.monomial( (j,i,ec) ) + if (j,i,ec) in A.indices(): + E_ij += A.monomial( (j,i,ec) ) else: - E_ij -= MS.monomial( (j,i,-ec) ) + E_ij -= A.monomial( (j,i,-ec) ) basis.append(E_ij) return tuple( basis ) @@ -2774,7 +2609,7 @@ class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): -2 """ - return (X*Y).trace().real().coefficient(0) + return (X*Y).trace().coefficient(0) class AlbertEJA(OctonionHermitianEJA): diff --git a/mjo/eja/eja_utils.py b/mjo/eja/eja_utils.py index 0b2d2a3..6f8cab6 100644 --- a/mjo/eja/eja_utils.py +++ b/mjo/eja/eja_utils.py @@ -54,7 +54,9 @@ def _all2list(x): SETUP:: sage: from mjo.eja.eja_utils import _all2list - sage: from mjo.octonions import Octonions, OctonionMatrixAlgebra + sage: from mjo.hurwitz import (QuaternionMatrixAlgebra, + ....: Octonions, + ....: OctonionMatrixAlgebra) EXAMPLES:: @@ -86,6 +88,13 @@ def _all2list(x): sage: _all2list(OctonionMatrixAlgebra(1).one()) [1, 0, 0, 0, 0, 0, 0, 0] + :: + + sage: _all2list(QuaternionAlgebra(QQ, -1, -1).one()) + [1, 0, 0, 0] + sage: _all2list(QuaternionMatrixAlgebra(1).one()) + [1, 0, 0, 0] + :: sage: V1 = VectorSpace(QQ,2) @@ -97,6 +106,13 @@ def _all2list(x): [3, 4, 1, 0, 0, 0, 0, 0, 0, 0] """ + if hasattr(x, 'list') and hasattr(x, 'to_vector'): + # This avoids calling to_vector() on a matrix algebra with + # e.g. quaternions where the returned vector is of the wrong + # length (three instead of four) because the quaternions don't + # know how many generators they have. + return _all2list(x.list()) + if hasattr(x, 'to_vector'): # This works on matrices of e.g. octonions directly, without # first needing to convert them to a list of octonions and -- 2.43.2