X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=5bb4cee119e31cde1fd17accc7cc6e0bdee7f041;hb=fcdcc7a37d4b5ee83d5ff0bc98fe2b63c61a7f51;hp=a8e16ec94aaa628f26533f654f33468f6a3eace4;hpb=17c1cf361b330252acae5ba18edb3a4fdf8bf9bd;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index a8e16ec..5bb4cee 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -1,9 +1,53 @@ """ -Euclidean Jordan Algebras. These are formally-real Jordan Algebras; -specifically those where u^2 + v^2 = 0 implies that u = v = 0. They -are used in optimization, and have some additional nice methods beyond -what can be supported in a general Jordan Algebra. - +Representations and constructions for Euclidean Jordan algebras. + +A Euclidean Jordan algebra is a Jordan algebra that has some +additional properties: + + 1. It is finite-dimensional. + 2. Its scalar field is the real numbers. + 3a. An inner product is defined on it, and... + 3b. That inner product is compatible with the Jordan product + in the sense that ` = ` for all elements + `x,y,z` in the algebra. + +Every Euclidean Jordan algebra is formally-real: for any two elements +`x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y = +0`. Conversely, every finite-dimensional formally-real Jordan algebra +can be made into a Euclidean Jordan algebra with an appropriate choice +of inner-product. + +Formally-real Jordan algebras were originally studied as a framework +for quantum mechanics. Today, Euclidean Jordan algebras are crucial in +symmetric cone optimization, since every symmetric cone arises as the +cone of squares in some Euclidean Jordan algebra. + +It is known that every Euclidean Jordan algebra decomposes into an +orthogonal direct sum (essentially, a Cartesian product) of simple +algebras, and that moreover, up to Jordan-algebra isomorphism, there +are only five families of simple algebras. We provide constructions +for these simple algebras: + + * :class:`BilinearFormEJA` + * :class:`RealSymmetricEJA` + * :class:`ComplexHermitianEJA` + * :class:`QuaternionHermitianEJA` + +Missing from this list is the algebra of three-by-three octononion +Hermitian matrices, as there is (as of yet) no implementation of the +octonions in SageMath. In addition to these, we provide two other +example constructions, + + * :class:`HadamardEJA` + * :class:`TrivialEJA` + +The Jordan spin algebra is a bilinear form algebra where the bilinear +form is the identity. The Hadamard EJA is simply a Cartesian product +of one-dimensional spin algebras. And last but not least, the trivial +EJA is exactly what you think. Cartesian products of these are also +supported using the usual ``cartesian_product()`` function; as a +result, we support (up to isomorphism) all Euclidean Jordan algebras +that don't involve octonions. SETUP:: @@ -13,7 +57,6 @@ EXAMPLES:: sage: random_eja() Euclidean Jordan algebra of dimension... - """ from itertools import repeat @@ -41,24 +84,50 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): INPUT: - - basis -- a tuple of basis elements in "matrix form," which - must be the same form as the arguments to ``jordan_product`` - and ``inner_product``. In reality, "matrix form" can be either - vectors, matrices, or a Cartesian product (ordered tuple) - of vectors or matrices. All of these would ideally be vector - spaces in sage with no special-casing needed; but in reality - we turn vectors into column-matrices and Cartesian products - `(a,b)` into column matrices `(a,b)^{T}` after converting - `a` and `b` themselves. - - - jordan_product -- function of two elements (in matrix form) - that returns their jordan product in this algebra; this will - be applied to ``basis`` to compute a multiplication table for - the algebra. - - - inner_product -- function of two elements (in matrix form) that - returns their inner product. This will be applied to ``basis`` to - compute an inner-product table (basically a matrix) for this algebra. + - ``basis`` -- a tuple; a tuple of basis elements in "matrix + form," which must be the same form as the arguments to + ``jordan_product`` and ``inner_product``. In reality, "matrix + form" can be either vectors, matrices, or a Cartesian product + (ordered tuple) of vectors or matrices. All of these would + ideally be vector spaces in sage with no special-casing + needed; but in reality we turn vectors into column-matrices + and Cartesian products `(a,b)` into column matrices + `(a,b)^{T}` after converting `a` and `b` themselves. + + - ``jordan_product`` -- a function; afunction of two ``basis`` + elements (in matrix form) that returns their jordan product, + also in matrix form; this will be applied to ``basis`` to + compute a multiplication table for the algebra. + + - ``inner_product`` -- a function; a function of two ``basis`` + elements (in matrix form) that returns their inner + product. This will be applied to ``basis`` to compute an + inner-product table (basically a matrix) for this algebra. + + - ``field`` -- a subfield of the reals (default: ``AA``); the scalar + field for the algebra. + + - ``orthonormalize`` -- boolean (default: ``True``); whether or + not to orthonormalize the basis. Doing so is expensive and + generally rules out using the rationals as your ``field``, but + is required for spectral decompositions. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS: + + We should compute that an element subalgebra is associative even + if we circumvent the element method:: + + sage: set_random_seed() + sage: J = random_eja(field=QQ,orthonormalize=False) + sage: x = J.random_element() + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: basis = tuple(b.superalgebra_element() for b in A.basis()) + sage: J.subalgebra(basis, orthonormalize=False).is_associative() + True """ Element = FiniteDimensionalEJAElement @@ -69,23 +138,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): inner_product, field=AA, orthonormalize=True, - associative=False, + associative=None, cartesian_product=False, check_field=True, check_axioms=True, prefix='e'): - # Keep track of whether or not the matrix basis consists of - # tuples, since we need special cases for them damned near - # everywhere. This is INDEPENDENT of whether or not the - # algebra is a cartesian product, since a subalgebra of a - # cartesian product will have a basis of tuples, but will not - # in general itself be a cartesian product algebra. - self._matrix_basis_is_cartesian = False n = len(basis) - if n > 0: - if hasattr(basis[0], 'cartesian_factors'): - self._matrix_basis_is_cartesian = True if check_field: if not field.is_subring(RR): @@ -94,20 +153,10 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # we've specified a real embedding. raise ValueError("scalar field is not real") + from mjo.eja.eja_utils import _change_ring # If the basis given to us wasn't over the field that it's # supposed to be over, fix that. Or, you know, crash. - if not cartesian_product: - # The field for a cartesian product algebra comes from one - # of its factors and is the same for all factors, so - # there's no need to "reapply" it on product algebras. - if self._matrix_basis_is_cartesian: - # OK since if n == 0, the basis does not consist of tuples. - P = basis[0].parent() - basis = tuple( P(tuple(b_i.change_ring(field) for b_i in b)) - for b in basis ) - else: - basis = tuple( b.change_ring(field) for b in basis ) - + basis = tuple( _change_ring(b, field) for b in basis ) if check_axioms: # Check commutativity of the Jordan and inner-products. @@ -126,12 +175,28 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): category = MagmaticAlgebras(field).FiniteDimensional() - category = category.WithBasis().Unital() + category = category.WithBasis().Unital().Commutative() + + 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() + # method because we need to know the category before we + # initialize the algebra. + associative = all( jordan_product(jordan_product(bi,bj),bk) + == + jordan_product(bi,jordan_product(bj,bk)) + for bi in basis + for bj in basis + for bk in basis) + if associative: # Element subalgebras can take advantage of this. category = category.Associative() if cartesian_product: - category = category.CartesianProducts() + # Use join() here because otherwise we only get the + # "Cartesian product of..." and not the things themselves. + category = category.join([category, + category.CartesianProducts()]) # Call the superclass constructor so that we can use its from_vector() # method to build our multiplication table. @@ -279,8 +344,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: if n > 0: ....: i = ZZ.random_element(n) ....: j = ZZ.random_element(n) - ....: ei = J.gens()[i] - ....: ej = J.gens()[j] + ....: ei = J.monomial(i) + ....: ej = J.monomial(j) ....: ei_ej = J.product_on_basis(i,j) sage: ei*ej == ei_ej True @@ -370,6 +435,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): """ return "Associative" in self.category().axioms() + def _is_commutative(self): + r""" + Whether or not this algebra's multiplication table is commutative. + + This method should of course always return ``True``, unless + this algebra was constructed with ``check_axioms=False`` and + passed an invalid multiplication table. + """ + return all( x*y == y*x for x in self.gens() for y in self.gens() ) + def _is_jordanian(self): r""" Whether or not this algebra's multiplication table respects the @@ -377,16 +452,102 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): We only check one arrangement of `x` and `y`, so for a ``True`` result to be truly true, you should also check - :meth:`is_commutative`. This method should of course always + :meth:`_is_commutative`. This method should of course always return ``True``, unless this algebra was constructed with ``check_axioms=False`` and passed an invalid multiplication table. """ - return all( (self.gens()[i]**2)*(self.gens()[i]*self.gens()[j]) + return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j)) == - (self.gens()[i])*((self.gens()[i]**2)*self.gens()[j]) + (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j)) for i in range(self.dimension()) for j in range(self.dimension()) ) + def _jordan_product_is_associative(self): + r""" + Return whether or not this algebra's Jordan product is + associative; that is, whether or not `x*(y*z) = (x*y)*z` + for all `x,y,x`. + + This method should agree with :meth:`is_associative` unless + you lied about the value of the ``associative`` parameter + when you constructed the algebra. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (random_eja, + ....: RealSymmetricEJA, + ....: ComplexHermitianEJA, + ....: QuaternionHermitianEJA) + + EXAMPLES:: + + sage: J = RealSymmetricEJA(4, orthonormalize=False) + sage: J._jordan_product_is_associative() + False + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by() + sage: A._jordan_product_is_associative() + True + + :: + + sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) + sage: J._jordan_product_is_associative() + False + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: A._jordan_product_is_associative() + True + + :: + + sage: J = QuaternionHermitianEJA(2) + sage: J._jordan_product_is_associative() + False + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by() + sage: A._jordan_product_is_associative() + True + + TESTS: + + The values we've presupplied to the constructors agree with + the computation:: + + sage: set_random_seed() + sage: J = random_eja() + sage: J.is_associative() == J._jordan_product_is_associative() + True + + """ + R = self.base_ring() + + # Used to check whether or not something is zero. + epsilon = R.zero() + if not R.is_exact(): + # I don't know of any examples that make this magnitude + # necessary because I don't know how to make an + # associative algebra when the element subalgebra + # construction is unreliable (as it is over RDF; we can't + # find the degree of an element because we can't compute + # the rank of a matrix). But even multiplication of floats + # is non-associative, so *some* epsilon is needed... let's + # just take the one from _inner_product_is_associative? + epsilon = 1e-15 + + for i in range(self.dimension()): + for j in range(self.dimension()): + for k in range(self.dimension()): + x = self.monomial(i) + y = self.monomial(j) + z = self.monomial(k) + diff = (x*y)*z - x*(y*z) + + if diff.norm() > epsilon: + return False + + return True + def _inner_product_is_associative(self): r""" Return whether or not this algebra's inner product `B` is @@ -396,26 +557,25 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): this algebra was constructed with ``check_axioms=False`` and passed an invalid Jordan or inner-product. """ + R = self.base_ring() - # Used to check whether or not something is zero in an inexact - # ring. This number is sufficient to allow the construction of - # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True. - epsilon = 1e-16 + # Used to check whether or not something is zero. + epsilon = R.zero() + if not R.is_exact(): + # This choice is sufficient to allow the construction of + # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True. + epsilon = 1e-15 for i in range(self.dimension()): for j in range(self.dimension()): for k in range(self.dimension()): - x = self.gens()[i] - y = self.gens()[j] - z = self.gens()[k] + x = self.monomial(i) + y = self.monomial(j) + z = self.monomial(k) diff = (x*y).inner_product(z) - x.inner_product(y*z) - if self.base_ring().is_exact(): - if diff != 0: - return False - else: - if diff.abs() > epsilon: - return False + if diff.abs() > epsilon: + return False return True @@ -429,7 +589,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + sage: from mjo.eja.eja_algebra import (random_eja, + ....: JordanSpinEJA, ....: HadamardEJA, ....: RealSymmetricEJA) @@ -458,22 +619,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: J2 = RealSymmetricEJA(2) sage: J = cartesian_product([J1,J2]) sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) ) - e(0, 1) + e(1, 2) + e1 + e5 TESTS: - Ensure that we can convert any element of the two non-matrix - simple algebras (whose matrix representations are columns) - back and forth faithfully:: + Ensure that we can convert any element back and forth + faithfully between its matrix and algebra representations:: sage: set_random_seed() - sage: J = HadamardEJA.random_instance() - sage: x = J.random_element() - sage: J(x.to_vector().column()) == x - True - sage: J = JordanSpinEJA.random_instance() + sage: J = random_eja() sage: x = J.random_element() - sage: J(x.to_vector().column()) == x + sage: J(x.to_matrix()) == x True We cannot coerce elements between algebras just because their @@ -489,7 +645,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Traceback (most recent call last): ... ValueError: not an element of this algebra - """ msg = "not an element of this algebra" if elt in self.base_ring(): @@ -757,7 +912,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # And to each subsequent row, prepend an entry that belongs to # the left-side "header column." - M += [ [self.gens()[i]] + [ self.product_on_basis(i,j) + M += [ [self.monomial(i)] + [ self.monomial(i)*self.monomial(j) for j in range(n) ] for i in range(n) ] @@ -828,12 +983,49 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): we think of them as matrices (including column vectors of the appropriate size). - Generally this will be an `n`-by-`1` column-vector space, + "By default" this will be an `n`-by-`1` column-matrix space, except when the algebra is trivial. There it's `n`-by-`n` (where `n` is zero), to ensure that two elements of the matrix - space (empty matrices) can be multiplied. + space (empty matrices) can be multiplied. For algebras of + matrices, this returns the space in which their + real embeddings live. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, + ....: JordanSpinEJA, + ....: QuaternionHermitianEJA, + ....: TrivialEJA) + + EXAMPLES: + + By default, the matrix representation is just a column-matrix + equivalent to the vector representation:: + + sage: J = JordanSpinEJA(3) + sage: J.matrix_space() + Full MatrixSpace of 3 by 1 dense matrices over Algebraic + Real Field + + The matrix representation in the trivial algebra is + zero-by-zero instead of the usual `n`-by-one:: + + sage: J = TrivialEJA() + sage: J.matrix_space() + Full MatrixSpace of 0 by 0 dense matrices over Algebraic + Real Field + + The matrix space for complex/quaternion Hermitian matrix EJA + is the space in which their real-embeddings live, not the + original complex/quaternion matrix space:: + + sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) + sage: J.matrix_space() + 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 - Matrix algebras override this with something more useful. """ if self.is_trivial(): return MatrixSpace(self.base_ring(), 0) @@ -1222,7 +1414,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): def L_x_i_j(i,j): # From a result in my book, these are the entries of the # basis representation of L_x. - return sum( vars[k]*self.gens()[k].operator().matrix()[i,j] + return sum( vars[k]*self.monomial(k).operator().matrix()[i,j] for k in range(n) ) L_x = matrix(F, n, n, L_x_i_j) @@ -1391,6 +1583,13 @@ class RationalBasisEJA(FiniteDimensionalEJA): if not all( all(b_i in QQ for b_i in b.list()) for b in basis ): raise TypeError("basis not rational") + super().__init__(basis, + jordan_product, + inner_product, + field=field, + check_field=check_field, + **kwargs) + self._rational_algebra = None if field is not QQ: # There's no point in constructing the extra algebra if this @@ -1404,17 +1603,11 @@ class RationalBasisEJA(FiniteDimensionalEJA): jordan_product, inner_product, field=QQ, + associative=self.is_associative(), orthonormalize=False, check_field=False, check_axioms=False) - super().__init__(basis, - jordan_product, - inner_product, - field=field, - check_field=check_field, - **kwargs) - @cached_method def _charpoly_coefficients(self): r""" @@ -1687,9 +1880,9 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): In theory, our "field" can be any subfield of the reals:: - sage: RealSymmetricEJA(2, field=RDF) + sage: RealSymmetricEJA(2, field=RDF, check_axioms=True) Euclidean Jordan algebra of dimension 3 over Real Double Field - sage: RealSymmetricEJA(2, field=RR) + sage: RealSymmetricEJA(2, field=RR, check_axioms=True) Euclidean Jordan algebra of dimension 3 over Real Field with 53 bits of precision @@ -1778,10 +1971,15 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): # if the user passes check_axioms=True. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - super(RealSymmetricEJA, self).__init__(self._denormalized_basis(n), - self.jordan_product, - self.trace_inner_product, - **kwargs) + associative = False + if n <= 1: + associative = True + + super().__init__(self._denormalized_basis(n), + self.jordan_product, + self.trace_inner_product, + associative=associative, + **kwargs) # TODO: this could be factored out somehow, but is left here # because the MatrixEJA is not presently a subclass of the @@ -1871,7 +2069,7 @@ class ComplexMatrixEJA(MatrixEJA): True """ - super(ComplexMatrixEJA,cls).real_embed(M) + super().real_embed(M) n = M.nrows() # We don't need any adjoined elements... @@ -1918,7 +2116,7 @@ class ComplexMatrixEJA(MatrixEJA): True """ - super(ComplexMatrixEJA,cls).real_unembed(M) + super().real_unembed(M) n = ZZ(M.nrows()) d = cls.dimension_over_reals() F = cls.complex_extension(M.base_ring()) @@ -1955,9 +2153,9 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): In theory, our "field" can be any subfield of the reals:: - sage: ComplexHermitianEJA(2, field=RDF) + sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True) Euclidean Jordan algebra of dimension 4 over Real Double Field - sage: ComplexHermitianEJA(2, field=RR) + sage: ComplexHermitianEJA(2, field=RR, check_axioms=True) Euclidean Jordan algebra of dimension 4 over Real Field with 53 bits of precision @@ -2066,10 +2264,15 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): # if the user passes check_axioms=True. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - super(ComplexHermitianEJA, self).__init__(self._denormalized_basis(n), - self.jordan_product, - self.trace_inner_product, - **kwargs) + associative = False + if n <= 1: + associative = True + + super().__init__(self._denormalized_basis(n), + self.jordan_product, + self.trace_inner_product, + 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(). @@ -2152,7 +2355,7 @@ class QuaternionMatrixEJA(MatrixEJA): True """ - super(QuaternionMatrixEJA,cls).real_embed(M) + super().real_embed(M) quaternions = M.base_ring() n = M.nrows() @@ -2207,7 +2410,7 @@ class QuaternionMatrixEJA(MatrixEJA): True """ - super(QuaternionMatrixEJA,cls).real_unembed(M) + super().real_unembed(M) n = ZZ(M.nrows()) d = cls.dimension_over_reals() @@ -2252,9 +2455,9 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA): In theory, our "field" can be any subfield of the reals:: - sage: QuaternionHermitianEJA(2, field=RDF) + sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True) Euclidean Jordan algebra of dimension 6 over Real Double Field - sage: QuaternionHermitianEJA(2, field=RR) + sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True) Euclidean Jordan algebra of dimension 6 over Real Field with 53 bits of precision @@ -2372,10 +2575,16 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA): # if the user passes check_axioms=True. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - super(QuaternionHermitianEJA, self).__init__(self._denormalized_basis(n), - self.jordan_product, - self.trace_inner_product, - **kwargs) + associative = False + if n <= 1: + associative = True + + super().__init__(self._denormalized_basis(n), + self.jordan_product, + self.trace_inner_product, + 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(). @@ -2601,10 +2810,17 @@ class BilinearFormEJA(ConcreteEJA): n = B.nrows() column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() ) - super(BilinearFormEJA, self).__init__(column_basis, - jordan_product, - inner_product, - **kwargs) + + # TODO: I haven't actually checked this, but it seems legit. + associative = False + if n <= 2: + associative = True + + super().__init__(column_basis, + jordan_product, + inner_product, + associative=associative, + **kwargs) # The rank of this algebra is two, unless we're in a # one-dimensional ambient space (because the rank is bounded @@ -2709,7 +2925,7 @@ class JordanSpinEJA(BilinearFormEJA): # But also don't pass check_field=False here, because the user # can pass in a field! - super(JordanSpinEJA, self).__init__(B, **kwargs) + super().__init__(B, **kwargs) @staticmethod def _max_random_instance_size(): @@ -2767,10 +2983,12 @@ class TrivialEJA(ConcreteEJA): if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - super(TrivialEJA, self).__init__(basis, - jordan_product, - inner_product, - **kwargs) + super().__init__(basis, + jordan_product, + inner_product, + associative=True, + **kwargs) + # The rank is zero using my definition, namely the dimension of the # largest subalgebra generated by any element. self.rank.set_cache(0) @@ -2783,8 +3001,7 @@ class TrivialEJA(ConcreteEJA): return cls(**kwargs) -class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct, - FiniteDimensionalEJA): +class CartesianProductEJA(FiniteDimensionalEJA): r""" The external (orthogonal) direct sum of two or more Euclidean Jordan algebras. Every Euclidean Jordan algebra decomposes into an @@ -2880,6 +3097,33 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct, sage: CP2.is_associative() False + Cartesian products of Cartesian products work:: + + sage: J1 = JordanSpinEJA(1) + sage: J2 = JordanSpinEJA(1) + sage: J3 = JordanSpinEJA(1) + sage: J = cartesian_product([J1,cartesian_product([J2,J3])]) + sage: J.multiplication_table() + +----++----+----+----+ + | * || e0 | e1 | e2 | + +====++====+====+====+ + | e0 || e0 | 0 | 0 | + +----++----+----+----+ + | e1 || 0 | e1 | 0 | + +----++----+----+----+ + | e2 || 0 | 0 | e2 | + +----++----+----+----+ + sage: HadamardEJA(3).multiplication_table() + +----++----+----+----+ + | * || e0 | e1 | e2 | + +====++====+====+====+ + | e0 || e0 | 0 | 0 | + +----++----+----+----+ + | e1 || 0 | e1 | 0 | + +----++----+----+----+ + | e2 || 0 | 0 | e2 | + +----++----+----+----+ + TESTS: All factors must share the same base field:: @@ -2907,37 +3151,41 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct, Element = FiniteDimensionalEJAElement - def __init__(self, algebras, **kwargs): - CombinatorialFreeModule_CartesianProduct.__init__(self, - algebras, - **kwargs) - field = algebras[0].base_ring() - if not all( J.base_ring() == field for J in algebras ): + def __init__(self, factors, **kwargs): + m = len(factors) + if m == 0: + return TrivialEJA() + + self._sets = factors + + field = factors[0].base_ring() + if not all( J.base_ring() == field for J in factors ): raise ValueError("all factors must share the same base field") - associative = all( m.is_associative() for m in algebras ) + associative = all( f.is_associative() for f in factors ) - # The definition of matrix_space() and self.basis() relies - # only on the stuff in the CFM_CartesianProduct class, which - # we've already initialized. - Js = self.cartesian_factors() - m = len(Js) MS = self.matrix_space() - basis = tuple( - MS(tuple( self.cartesian_projection(i)(b).to_matrix() - for i in range(m) )) - for b in self.basis() - ) + basis = [] + zero = MS.zero() + for i in range(m): + for b in factors[i].matrix_basis(): + z = list(zero) + z[i] = b + basis.append(z) + + basis = tuple( MS(b) for b in basis ) # Define jordan/inner products that operate on that matrix_basis. def jordan_product(x,y): return MS(tuple( - (Js[i](x[i])*Js[i](y[i])).to_matrix() for i in range(m) + (factors[i](x[i])*factors[i](y[i])).to_matrix() + for i in range(m) )) def inner_product(x, y): return sum( - Js[i](x[i]).inner_product(Js[i](y[i])) for i in range(m) + factors[i](x[i]).inner_product(factors[i](y[i])) + for i in range(m) ) # There's no need to check the field since it already came @@ -2957,9 +3205,25 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct, check_field=False, check_axioms=False) - ones = tuple(J.one() for J in algebras) - self.one.set_cache(self._cartesian_product_of_elements(ones)) - self.rank.set_cache(sum(J.rank() for J in algebras)) + ones = tuple(J.one().to_matrix() for J in factors) + self.one.set_cache(self(ones)) + self.rank.set_cache(sum(J.rank() for J in factors)) + + def cartesian_factors(self): + # Copy/pasted from CombinatorialFreeModule_CartesianProduct. + return self._sets + + def cartesian_factor(self, i): + r""" + Return the ``i``th factor of this algebra. + """ + return self._sets[i] + + def _repr_(self): + # Copy/pasted from CombinatorialFreeModule_CartesianProduct. + from sage.categories.cartesian_product import cartesian_product + return cartesian_product.symbol.join("%s" % factor + for factor in self._sets) def matrix_space(self): r""" @@ -3056,9 +3320,12 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct, True """ - Ji = self.cartesian_factors()[i] - # Requires the fix on Trac 31421/31422 to work! - Pi = super().cartesian_projection(i) + offset = sum( self.cartesian_factor(k).dimension() + for k in range(i) ) + Ji = self.cartesian_factor(i) + Pi = self._module_morphism(lambda j: Ji.monomial(j - offset), + codomain=Ji) + return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix()) @cached_method @@ -3164,13 +3431,64 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct, True """ - Ji = self.cartesian_factors()[i] - # Requires the fix on Trac 31421/31422 to work! - Ei = super().cartesian_embedding(i) + offset = sum( self.cartesian_factor(k).dimension() + for k in range(i) ) + Ji = self.cartesian_factor(i) + Ei = Ji._module_morphism(lambda j: self.monomial(j + offset), + codomain=self) return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix()) FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA -random_eja = ConcreteEJA.random_instance +class RationalBasisCartesianProductEJA(CartesianProductEJA, + RationalBasisEJA): + r""" + A separate class for products of algebras for which we know a + rational basis. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + ....: RealSymmetricEJA) + + EXAMPLES: + + This gives us fast characteristic polynomial computations in + product algebras, too:: + + + sage: J1 = JordanSpinEJA(2) + sage: J2 = RealSymmetricEJA(3) + sage: J = cartesian_product([J1,J2]) + sage: J.characteristic_polynomial_of().degree() + 5 + sage: J.rank() + 5 + + """ + def __init__(self, algebras, **kwargs): + CartesianProductEJA.__init__(self, algebras, **kwargs) + + self._rational_algebra = None + if self.vector_space().base_field() is not QQ: + self._rational_algebra = cartesian_product([ + r._rational_algebra for r in algebras + ]) + + +RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA + +def random_eja(*args, **kwargs): + J1 = ConcreteEJA.random_instance(*args, **kwargs) + + # This might make Cartesian products appear roughly as often as + # any other ConcreteEJA. + if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0: + # Use random_eja() again so we can get more than two factors. + J2 = random_eja(*args, **kwargs) + J = cartesian_product([J1,J2]) + return J + else: + return J1