From 251e308af6c3a57b418a06546a572af83eade7ec Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Tue, 1 Dec 2020 10:54:27 -0500 Subject: [PATCH] eja: convert matrix algebas to the new constructor, fix all tests. --- mjo/eja/TODO | 7 + mjo/eja/eja_algebra.py | 264 ++++++++++-------------------- mjo/eja/eja_element.py | 37 ++++- mjo/eja/eja_element_subalgebra.py | 4 +- 4 files changed, 133 insertions(+), 179 deletions(-) diff --git a/mjo/eja/TODO b/mjo/eja/TODO index 13cd7f9..cbe1d95 100644 --- a/mjo/eja/TODO +++ b/mjo/eja/TODO @@ -48,3 +48,10 @@ 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 12. Don't pass in an n-by-n multiplication/i-p table since only the lower-left half is used. + +13. Move the "field" argument to a keyword after basis, jp, and ip. + +14. Instead of storing the multiplication and inner-product tables in + RationalBasisEuclideanJordanAlgebra, why not just create the + algebra over QQ in the constructor and save that? They're globally + unique, so we won't wind up with multiple copies. diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index effafd4..efa10e2 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -478,7 +478,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = HadamardEJA(2) sage: J.coordinate_polynomial_ring() Multivariate Polynomial Ring in X1, X2... - sage: J = RealSymmetricEJA(3,QQ) + sage: J = RealSymmetricEJA(3,QQ,orthonormalize=False) sage: J.coordinate_polynomial_ring() Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6... @@ -1163,6 +1163,7 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr if n > 0: if is_Matrix(basis[0]): basis_is_matrices = True + from mjo.eja.eja_utils import _vec2mat vector_basis = tuple( map(_mat2vec,basis) ) degree = basis[0].nrows()**2 else: @@ -1223,25 +1224,43 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr if self._deortho_inner_product_table is not None: self._deortho_inner_product_table = tuple(map(tuple, self._deortho_inner_product_table)) + # We overwrite the name "vector_basis" in a second, but never modify it + # in place, to this effectively makes a copy of it. + deortho_vector_basis = vector_basis + self._deortho_matrix = None + if orthonormalize: from mjo.eja.eja_utils import gram_schmidt - vector_basis = gram_schmidt(vector_basis, inner_product) + if basis_is_matrices: + vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y)) + vector_basis = gram_schmidt(vector_basis, vector_ip) + else: + vector_basis = gram_schmidt(vector_basis, inner_product) + W = V.span_of_basis( vector_basis ) # Normalize the "matrix" basis, too! basis = vector_basis if basis_is_matrices: - from mjo.eja.eja_utils import _vec2mat basis = tuple( map(_vec2mat,basis) ) W = V.span_of_basis( vector_basis ) + # Now "W" is the vector space of our algebra coordinates. The + # variables "X1", "X2",... refer to the entries of vectors in + # W. Thus to convert back and forth between the orthonormal + # coordinates and the given ones, we need to stick the original + # basis in W. + U = V.span_of_basis( deortho_vector_basis ) + self._deortho_matrix = matrix( U.coordinate_vector(q) + for q in vector_basis ) + # TODO: use symmetry mult_table = [ [0 for j in range(n)] for i in range(n) ] ip_table = [ [0 for j in range(n)] for i in range(n) ] - # Note: the Jordan and inner- products are defined in terms + # Note: the Jordan and inner-products are defined in terms # of the ambient basis. It's important that their arguments # are in ambient coordinates as well. for i in range(n): @@ -1294,19 +1313,6 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr EXAMPLES: - The returned coefficients should be the same as if we'd never - orthonormalized the basis to begin with:: - - sage: B = matrix(QQ, [[1, 0, 0], - ....: [0, 25, -32], - ....: [0, -32, 41] ]) - sage: J1 = BilinearFormEJA(B) - sage: J2 = BilinearFormEJA(B,QQ,orthonormalize=False) - sage: J1._charpoly_coefficients() - (X1^2 - 25*X2^2 + 64*X2*X3 - 41*X3^2, -2*X1) - sage: J2._charpoly_coefficients() - (X1^2 - 25*X2^2 + 64*X2*X3 - 41*X3^2, -2*X1) - The base ring of the resulting polynomial coefficients is what it should be, and not the rationals (unless the algebra was already over the rationals):: @@ -1319,6 +1325,7 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr Algebraic Real Field sage: a0.base_ring() Algebraic Real Field + """ if self.base_ring() is QQ: # There's no need to construct *another* algebra over the @@ -1331,10 +1338,22 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr J = FiniteDimensionalEuclideanJordanAlgebra(QQ, self._deortho_multiplication_table, self._deortho_inner_product_table) - a = J._charpoly_coefficients() - return tuple(map(lambda x: x.change_ring(self.base_ring()), a)) -class ConcreteEuclideanJordanAlgebra: + # Change back from QQ to our real base ring + a = ( a_i.change_ring(self.base_ring()) + for a_i in J._charpoly_coefficients() ) + + # Now convert the coordinate variables back to the + # deorthonormalized ones. + R = self.coordinate_polynomial_ring() + from sage.modules.free_module_element import vector + X = vector(R, R.gens()) + BX = self._deortho_matrix*X + + subs_dict = { X[i]: BX[i] for i in range(len(X)) } + return tuple( a_i.subs(subs_dict) for a_i in a ) + +class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra): r""" A class for the Euclidean Jordan algebras that we know by name. @@ -1387,7 +1406,7 @@ class ConcreteEuclideanJordanAlgebra: raise NotImplementedError @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, *args, **kwargs): """ Return a random instance of this type of algebra. @@ -1395,127 +1414,14 @@ class ConcreteEuclideanJordanAlgebra: """ from sage.misc.prandom import choice eja_class = choice(cls.__subclasses__()) - return eja_class.random_instance(field) - - -class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): - - def __init__(self, field, basis, normalize_basis=True, **kwargs): - """ - Compared to the superclass constructor, we take a basis instead of - a multiplication table because the latter can be computed in terms - of the former when the product is known (like it is here). - """ - # Used in this class's fast _charpoly_coefficients() override. - self._basis_normalizers = None - - # We're going to loop through this a few times, so now's a good - # time to ensure that it isn't a generator expression. - basis = tuple(basis) - - algebra_dim = len(basis) - degree = 0 # size of the matrices - if algebra_dim > 0: - degree = basis[0].nrows() - - if algebra_dim > 1 and normalize_basis: - # We'll need sqrt(2) to normalize the basis, and this - # winds up in the multiplication table, so the whole - # algebra needs to be over the field extension. - R = PolynomialRing(field, 'z') - z = R.gen() - p = z**2 - 2 - if p.is_irreducible(): - field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt()) - basis = tuple( s.change_ring(field) for s in basis ) - self._basis_normalizers = tuple( - ~(self.matrix_inner_product(s,s).sqrt()) for s in basis ) - basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers)) - - # Now compute the multiplication and inner product tables. - # We have to do this *after* normalizing the basis, because - # scaling affects the answers. - V = VectorSpace(field, degree**2) - W = V.span_of_basis( _mat2vec(s) for s in basis ) - mult_table = [[W.zero() for j in range(algebra_dim)] - for i in range(algebra_dim)] - ip_table = [[field.zero() for j in range(algebra_dim)] - for i in range(algebra_dim)] - for i in range(algebra_dim): - for j in range(algebra_dim): - mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2 - mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry)) - - try: - # HACK: ignore the error here if we don't need the - # inner product (as is the case when we construct - # a dummy QQ-algebra for fast charpoly coefficients. - ip_table[i][j] = self.matrix_inner_product(basis[i], - basis[j]) - except: - pass - - super(MatrixEuclideanJordanAlgebra, self).__init__(field, - mult_table, - ip_table, - matrix_basis=basis, - **kwargs) - - if algebra_dim == 0: - self.one.set_cache(self.zero()) - else: - n = basis[0].nrows() - # The identity wrt (A,B) -> (AB + BA)/2 is independent of the - # details of this algebra. - self.one.set_cache(self(matrix.identity(field,n))) - - @cached_method - def _charpoly_coefficients(self): - r""" - Override the parent method with something that tries to compute - over a faster (non-extension) field. - """ - if self._basis_normalizers is None or self.base_ring() is QQ: - # We didn't normalize, or the basis we started with had - # entries in a nice field already. Just compute the thing. - return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients() - - basis = ( (b/n) for (b,n) in zip(self.matrix_basis(), - self._basis_normalizers) ) - - # Do this over the rationals and convert back at the end. - # Only works because we know the entries of the basis are - # integers. The argument ``check_axioms=False`` is required - # because the trace inner-product method for this - # class is a stub and can't actually be checked. - J = MatrixEuclideanJordanAlgebra(QQ, - basis, - normalize_basis=False, - check_field=False, - check_axioms=False) - a = J._charpoly_coefficients() - - # Unfortunately, changing the basis does change the - # coefficients of the characteristic polynomial, but since - # these are really the coefficients of the "characteristic - # polynomial of" function, everything is still nice and - # unevaluated. It's therefore "obvious" how scaling the - # basis affects the coordinate variables X1, X2, et - # cetera. Scaling the first basis vector up by "n" adds a - # factor of 1/n into every "X1" term, for example. So here - # we simply undo the basis_normalizer scaling that we - # performed earlier. - # - # The a[0] access here is safe because trivial algebras - # won't have any basis normalizers and therefore won't - # make it to this "else" branch. - XS = a[0].parent().gens() - subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i] - for i in range(len(XS)) } - return tuple( a_i.subs(subs_dict) for a_i in a ) + # These all bubble up to the RationalBasisEuclideanJordanAlgebra + # superclass constructor, so any (kw)args valid there are also + # valid here. + return eja_class.random_instance(*args, **kwargs) +class MatrixEuclideanJordanAlgebra: @staticmethod def real_embed(M): """ @@ -1540,8 +1446,12 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): """ raise NotImplementedError + @staticmethod + def jordan_product(X,Y): + return (X*Y + Y*X)/2 + @classmethod - def matrix_inner_product(cls,X,Y): + def trace_inner_product(cls,X,Y): Xu = cls.real_unembed(X) Yu = cls.real_unembed(Y) tr = (Xu*Yu).trace() @@ -1574,8 +1484,8 @@ class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): return M -class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, - ConcreteEuclideanJordanAlgebra): +class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra, + RealMatrixEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -1688,9 +1598,11 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, basis = self._denormalized_basis(n, field) super(RealSymmetricEJA, self).__init__(field, basis, - check_axioms=False, + self.jordan_product, + self.trace_inner_product, **kwargs) self.rank.set_cache(n) + self.one.set_cache(self(matrix.identity(field,n))) class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @@ -1821,7 +1733,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @classmethod - def matrix_inner_product(cls,X,Y): + def trace_inner_product(cls,X,Y): """ Compute a matrix inner product in this algebra directly from its real embedding. @@ -1843,16 +1755,16 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): sage: X = ComplexHermitianEJA.real_unembed(Xe) sage: Y = ComplexHermitianEJA.real_unembed(Ye) sage: expected = (X*Y).trace().real() - sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye) + sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye) sage: actual == expected True """ - return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/2 + return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2 -class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, - ConcreteEuclideanJordanAlgebra): +class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, + ComplexMatrixEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, @@ -1961,16 +1873,18 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, # Since we embedded these, we can drop back to the "field" that we # started with instead of the complex extension "F". - return ( s.change_ring(field) for s in S ) + return tuple( s.change_ring(field) for s in S ) def __init__(self, n, field=AA, **kwargs): basis = self._denormalized_basis(n,field) - super(ComplexHermitianEJA,self).__init__(field, - basis, - check_axioms=False, - **kwargs) + super(ComplexHermitianEJA, self).__init__(field, + basis, + self.jordan_product, + self.trace_inner_product, + **kwargs) self.rank.set_cache(n) + # TODO: pre-cache the identity! @staticmethod def _max_random_instance_size(): @@ -2116,7 +2030,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @classmethod - def matrix_inner_product(cls,X,Y): + def trace_inner_product(cls,X,Y): """ Compute a matrix inner product in this algebra directly from its real embedding. @@ -2138,16 +2052,16 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): sage: X = QuaternionHermitianEJA.real_unembed(Xe) sage: Y = QuaternionHermitianEJA.real_unembed(Ye) sage: expected = (X*Y).trace().coefficient_tuple()[0] - sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye) + sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye) sage: actual == expected True """ - return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/4 + return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4 -class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, - ConcreteEuclideanJordanAlgebra): +class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, + QuaternionMatrixEuclideanJordanAlgebra): r""" The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the @@ -2257,16 +2171,18 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, # Since we embedded these, we can drop back to the "field" that we # started with instead of the quaternion algebra "Q". - return ( s.change_ring(field) for s in S ) + return tuple( s.change_ring(field) for s in S ) def __init__(self, n, field=AA, **kwargs): basis = self._denormalized_basis(n,field) - super(QuaternionHermitianEJA,self).__init__(field, - basis, - check_axioms=False, - **kwargs) + super(QuaternionHermitianEJA, self).__init__(field, + basis, + self.jordan_product, + self.trace_inner_product, + **kwargs) self.rank.set_cache(n) + # TODO: cache one()! @staticmethod def _max_random_instance_size(): @@ -2284,8 +2200,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, return cls(n, field, **kwargs) -class HadamardEJA(RationalBasisEuclideanJordanAlgebra, - ConcreteEuclideanJordanAlgebra): +class HadamardEJA(ConcreteEuclideanJordanAlgebra): """ Return the Euclidean Jordan Algebra corresponding to the set `R^n` under the Hadamard product. @@ -2362,8 +2277,7 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra, return cls(n, field, **kwargs) -class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra, - ConcreteEuclideanJordanAlgebra): +class BilinearFormEJA(ConcreteEuclideanJordanAlgebra): r""" The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)`` with the half-trace inner product and jordan product ``x*y = @@ -2585,8 +2499,7 @@ class JordanSpinEJA(BilinearFormEJA): return cls(n, field, **kwargs) -class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, - ConcreteEuclideanJordanAlgebra): +class TrivialEJA(ConcreteEuclideanJordanAlgebra): """ The trivial Euclidean Jordan algebra consisting of only a zero element. @@ -2616,12 +2529,13 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, """ def __init__(self, field=AA, **kwargs): - mult_table = [] - ip_table = [] + jordan_product = lambda x,y: x + inner_product = lambda x,y: field(0) + basis = () super(TrivialEJA, self).__init__(field, - mult_table, - ip_table, - check_axioms=False, + basis, + jordan_product, + inner_product, **kwargs) # The rank is zero using my definition, namely the dimension of the # largest subalgebra generated by any element. @@ -2849,7 +2763,7 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): EXAMPLE:: sage: J1 = HadamardEJA(3,QQ) - sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False) + sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=False) sage: J = DirectSumEJA(J1,J2) sage: x1 = J1.one() sage: x2 = x1 diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py index 6547668..0fae507 100644 --- a/mjo/eja/eja_element.py +++ b/mjo/eja/eja_element.py @@ -340,6 +340,8 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: TrivialEJA, + ....: RealSymmetricEJA, + ....: ComplexHermitianEJA, ....: random_eja) EXAMPLES:: @@ -387,6 +389,37 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: x,y = J.random_elements(2) sage: (x*y).det() == x.det()*y.det() True + + The determinant in matrix algebras is just the usual determinant:: + + sage: set_random_seed() + sage: X = matrix.random(QQ,3) + sage: X = X + X.T + sage: J1 = RealSymmetricEJA(3) + sage: J2 = RealSymmetricEJA(3,QQ,orthonormalize=False) + sage: expected = X.det() + sage: actual1 = J1(X).det() + sage: actual2 = J2(X).det() + sage: actual1 == expected + True + sage: actual2 == expected + True + + :: + + sage: set_random_seed() + sage: J1 = ComplexHermitianEJA(3) + sage: J2 = ComplexHermitianEJA(3,field=QQ,orthonormalize=False) + sage: X = matrix.random(GaussianIntegers(),3) + sage: X = X + X.H + sage: expected = AA(X.det()) + sage: actual1 = J1(J1.real_embed(X)).det() + sage: actual2 = J2(J2.real_embed(X)).det() + sage: expected == actual1 + True + sage: expected == actual2 + True + """ P = self.parent() r = P.rank() @@ -523,7 +556,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): zero, but we need the characteristic polynomial for the determinant. The minimal polynomial is a lot easier to get, so we use Corollary 2 in Chapter V of Koecher to check - whether or not the paren't algebra's zero element is a root + whether or not the parent algebra's zero element is a root of this element's minimal polynomial. That is... unless the coefficients of our algebra's @@ -945,7 +978,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: n_max = RealSymmetricEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) sage: J1 = RealSymmetricEJA(n) - sage: J2 = RealSymmetricEJA(n,normalize_basis=False) + sage: J2 = RealSymmetricEJA(n,orthonormalize=False) sage: X = random_matrix(AA,n) sage: X = X*X.transpose() sage: x1 = J1(X) diff --git a/mjo/eja/eja_element_subalgebra.py b/mjo/eja/eja_element_subalgebra.py index 7bc4b3a..7edf1df 100644 --- a/mjo/eja/eja_element_subalgebra.py +++ b/mjo/eja/eja_element_subalgebra.py @@ -105,7 +105,7 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide The identity element acts like the identity over the rationals:: sage: set_random_seed() - sage: x = random_eja(field=QQ).random_element() + sage: x = random_eja(field=QQ,orthonormalize=False).random_element() sage: A = x.subalgebra_generated_by() sage: x = A.random_element() sage: A.one()*x == x and x*A.one() == x @@ -125,7 +125,7 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide the rationals:: sage: set_random_seed() - sage: x = random_eja(field=QQ).random_element() + sage: x = random_eja(field=QQ,orthonormalize=False).random_element() sage: A = x.subalgebra_generated_by() sage: actual = A.one().operator().matrix() sage: expected = matrix.identity(A.base_ring(), A.dimension()) -- 2.44.2