From: Michael Orlitzky Date: Thu, 3 Dec 2020 19:24:44 +0000 (-0500) Subject: eja: refactor some matrix algebra stuff and break the tests. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=ea99a21239882d478c1a458a0411b1eb0588b84b;p=sage.d.git eja: refactor some matrix algebra stuff and break the tests. --- diff --git a/mjo/eja/TODO b/mjo/eja/TODO index 427a953..cf870ae 100644 --- a/mjo/eja/TODO +++ b/mjo/eja/TODO @@ -15,10 +15,5 @@ RealSymmetricEJA(4): sage: F = J.base_ring() sage: a0 = (1/4)*X[4]**2*X[6]**2 - (1/2)*X[2]*X[5]*X[6]**2 - (1/2)*X[3]*X[4]*X[6]*X[7] + (F(2).sqrt()/2)*X[1]*X[5]*X[6]*X[7] + (1/4)*X[3]**2*X[7]**2 - (1/2)*X[0]*X[5]*X[7]**2 + (F(2).sqrt()/2)*X[2]*X[3]*X[6]*X[8] - (1/2)*X[1]*X[4]*X[6*X[8] - (1/2)*X[1]*X[3]*X[7]*X[8] + (F(2).sqrt()/2)*X[0]*X[4]*X[7]*X[8] + (1/4)*X[1]**2*X[8]**2 - (1/2)*X[0]*X[2]*X[8]**2 - (1/2)*X[2]*X[3]**2*X[9] + (F(2).sqrt()/2)*X[1]*X[3]*X[4]*X[9] - (1/2)*X[0]*X[4]**2*X[9] - (1/2)*X[1]**2*X[5]*X[9] + X[0]*X[2]*X[5]*X[9] -5. Compute the scalar in the general natural_inner_product() for - matrices, so no overrides are necessary. Actually, this is - probably better implemented as a dimension_over_reals() method - that returns 1, 2, or 4. - -6. The main EJA element constructor is happy to convert between +5. The main EJA element constructor is happy to convert between e.g. HadamardEJA(3) and JordanSpinEJA(3). diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index e075ed2..5dee5d0 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -1459,7 +1459,22 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra): class MatrixEuclideanJordanAlgebra: @staticmethod - def real_embed(M): + def dimension_over_reals(): + r""" + The dimension of this matrix's base ring over the reals. + + 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. + + This is used to determine the size of the matrix returned from + :meth:`real_embed`, among other things. + """ + raise NotImplementedError + + @classmethod + def real_embed(cls,M): """ Embed the matrix ``M`` into a space of real matrices. @@ -1472,15 +1487,21 @@ class MatrixEuclideanJordanAlgebra: real_embed(M*N) = real_embed(M)*real_embed(N) """ - raise NotImplementedError + if M.ncols() != M.nrows(): + raise ValueError("the matrix 'M' must be square") + return M - @staticmethod - def real_unembed(M): + @classmethod + def real_unembed(cls,M): """ The inverse of :meth:`real_embed`. """ - raise NotImplementedError + 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 @staticmethod def jordan_product(X,Y): @@ -1488,36 +1509,65 @@ class MatrixEuclideanJordanAlgebra: @classmethod def trace_inner_product(cls,X,Y): + r""" + Compute the trace inner-product of two real-embeddings. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (RealSymmetricEJA, + ....: ComplexHermitianEJA, + ....: QuaternionHermitianEJA) + + EXAMPLES:: + + This gives the same answer as it would if we computed the trace + from the unembedded (original) matrices:: + + 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 = ComplexHermitianEJA.real_unembed(Xe) + sage: Y = ComplexHermitianEJA.real_unembed(Ye) + sage: expected = (X*Y).trace().real() + sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye) + 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 = QuaternionHermitianEJA.real_unembed(Xe) + sage: Y = QuaternionHermitianEJA.real_unembed(Ye) + sage: expected = (X*Y).trace().coefficient_tuple()[0] + sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye) + sage: actual == expected + True + + """ Xu = cls.real_unembed(X) Yu = cls.real_unembed(Y) tr = (Xu*Yu).trace() try: # Works in QQ, AA, RDF, et cetera. - return tr.real() + return tr.real() / cls.dimension_over_reals() except AttributeError: # A quaternion doesn't have a real() method, but does # have coefficient_tuple() method that returns the # coefficients of 1, i, j, and k -- in that order. - return tr.coefficient_tuple()[0] + return tr.coefficient_tuple()[0] / cls.dimension_over_reals() class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @staticmethod - def real_embed(M): - """ - The identity function, for embedding real matrices into real - matrices. - """ - return M - - @staticmethod - def real_unembed(M): - """ - The identity function, for unembedding real matrices from real - matrices. - """ - return M + def dimension_over_reals(): + return 1 class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra, @@ -1639,13 +1689,23 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra, self.jordan_product, self.trace_inner_product, **kwargs) + + # TODO: this could be factored out somehow, but is left here + # because the MatrixEuclideanJordanAlgebra is not presently + # a subclass of the FDEJA class that defines rank() and one(). self.rank.set_cache(n) - self.one.set_cache(self(matrix.identity(ZZ,n))) + idV = matrix.identity(ZZ, self.dimension_over_reals()*n) + self.one.set_cache(self(idV)) + class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @staticmethod - def real_embed(M): + 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 + @@ -1687,9 +1747,8 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): True """ + super(ComplexMatrixEuclideanJordanAlgebra,cls).real_embed(M) n = M.nrows() - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") # We don't need any adjoined elements... field = M.base_ring().base_ring() @@ -1703,8 +1762,8 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): return matrix.block(field, n, blocks) - @staticmethod - def real_unembed(M): + @classmethod + def real_unembed(cls,M): """ The inverse of _embed_complex_matrix(). @@ -1735,11 +1794,9 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): True """ + super(ComplexMatrixEuclideanJordanAlgebra,cls).real_unembed(M) n = ZZ(M.nrows()) - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - if not n.mod(2).is_zero(): - raise ValueError("the matrix 'M' must be a complex embedding") + d = cls.dimension_over_reals() # If "M" was normalized, its base ring might have roots # adjoined and they can stick around after unembedding. @@ -1757,9 +1814,9 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): # 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/2): - for j in range(n/2): - submat = M[2*k:2*k+2,2*j:2*j+2] + 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]: @@ -1767,38 +1824,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): z = submat[0,0] + submat[0,1]*i elements.append(z) - return matrix(F, n/2, elements) - - - @classmethod - def trace_inner_product(cls,X,Y): - """ - Compute a matrix 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.to_matrix() - sage: Ye = y.to_matrix() - sage: X = ComplexHermitianEJA.real_unembed(Xe) - sage: Y = ComplexHermitianEJA.real_unembed(Ye) - sage: expected = (X*Y).trace().real() - sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye) - sage: actual == expected - True - - """ - return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2 + return matrix(F, n/d, elements) class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, @@ -1924,8 +1950,12 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, self.jordan_product, self.trace_inner_product, **kwargs) + # TODO: this could be factored out somehow, but is left here + # because the MatrixEuclideanJordanAlgebra is not presently + # a subclass of the FDEJA class that defines rank() and one(). self.rank.set_cache(n) - # TODO: pre-cache the identity! + idV = matrix.identity(ZZ, self.dimension_over_reals()*n) + self.one.set_cache(self(idV)) @staticmethod def _max_random_instance_size(): @@ -1941,7 +1971,11 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @staticmethod - def real_embed(M): + 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 @@ -1980,10 +2014,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): True """ + super(QuaternionMatrixEuclideanJordanAlgebra,cls).real_embed(M) quaternions = M.base_ring() n = M.nrows() - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") F = QuadraticField(-1, 'I') i = F.gen() @@ -2006,8 +2039,8 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): - @staticmethod - def real_unembed(M): + @classmethod + def real_unembed(cls,M): """ The inverse of _embed_quaternion_matrix(). @@ -2037,11 +2070,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): True """ + super(QuaternionMatrixEuclideanJordanAlgebra,cls).real_unembed(M) n = ZZ(M.nrows()) - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - if not n.mod(4).is_zero(): - raise ValueError("the matrix 'M' must be a quaternion embedding") + 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. @@ -2053,10 +2084,10 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): # 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/4): - for m in range(n/4): + for l in range(n/d): + for m in range(n/d): submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed( - M[4*l:4*l+4,4*m:4*m+4] ) + 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(): @@ -2067,38 +2098,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): z += submat[0,1].imag()*k elements.append(z) - return matrix(Q, n/4, elements) - - - @classmethod - def trace_inner_product(cls,X,Y): - """ - Compute a matrix 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.to_matrix() - sage: Ye = y.to_matrix() - sage: X = QuaternionHermitianEJA.real_unembed(Xe) - sage: Y = QuaternionHermitianEJA.real_unembed(Ye) - sage: expected = (X*Y).trace().coefficient_tuple()[0] - sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye) - sage: actual == expected - True - - """ - return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4 + return matrix(Q, n/d, elements) class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, @@ -2225,8 +2225,13 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, self.jordan_product, self.trace_inner_product, **kwargs) + # TODO: this could be factored out somehow, but is left here + # because the MatrixEuclideanJordanAlgebra is not presently + # a subclass of the FDEJA class that defines rank() and one(). self.rank.set_cache(n) - # TODO: cache one()! + idV = matrix.identity(ZZ, self.dimension_over_reals()*n) + self.one.set_cache(self(idV)) + @staticmethod def _max_random_instance_size():