X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=48421e357544ef2b5c5613bf9e222827e9a20fec;hb=5b1f16399286eba471884a1cfe45247b3a0a7693;hp=7c8adc7400880b0c4a95f7d73dec6cd11113aec6;hpb=a4f0908c2216ff989161d33873102805d1c6aabd;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 7c8adc7..48421e3 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,14 +57,15 @@ EXAMPLES:: sage: random_eja() Euclidean Jordan algebra of dimension... - """ from itertools import repeat from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra from sage.categories.magmatic_algebras import MagmaticAlgebras -from sage.combinat.free_module import CombinatorialFreeModule +from sage.categories.sets_cat import cartesian_product +from sage.combinat.free_module import (CombinatorialFreeModule, + CombinatorialFreeModule_CartesianProduct) from sage.matrix.constructor import matrix from sage.matrix.matrix_space import MatrixSpace from sage.misc.cachefunc import cached_method @@ -31,7 +76,7 @@ from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF, QuadraticField) from mjo.eja.eja_element import FiniteDimensionalEJAElement from mjo.eja.eja_operator import FiniteDimensionalEJAOperator -from mjo.eja.eja_utils import _mat2vec +from mjo.eja.eja_utils import _all2list, _mat2vec class FiniteDimensionalEJA(CombinatorialFreeModule): r""" @@ -39,16 +84,33 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): INPUT: - - basis -- a tuple of basis elements in their matrix form. - - - 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. """ Element = FiniteDimensionalEJAElement @@ -60,10 +122,23 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): field=AA, orthonormalize=True, associative=False, + 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): # Note: this does return true for the real algebraic @@ -73,7 +148,18 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # If the basis given to us wasn't over the field that it's # supposed to be over, fix that. Or, you know, crash. - basis = tuple( b.change_ring(field) for b in basis ) + 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 ) + if check_axioms: # Check commutativity of the Jordan and inner-products. @@ -96,15 +182,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): if associative: # Element subalgebras can take advantage of this. category = category.Associative() + if cartesian_product: + category = category.CartesianProducts() # Call the superclass constructor so that we can use its from_vector() # method to build our multiplication table. - n = len(basis) - super().__init__(field, - range(n), - prefix=prefix, - category=category, - bracket=False) + CombinatorialFreeModule.__init__(self, + field, + range(n), + prefix=prefix, + category=category, + bracket=False) # Now comes all of the hard work. We'll be constructing an # ambient vector space V that our (vectorized) basis lives in, @@ -115,8 +203,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): degree = 0 if n > 0: - # Works on both column and square matrices... - degree = len(basis[0].list()) + degree = len(_all2list(basis[0])) # Build an ambient space that fits our matrix basis when # written out as "long vectors." @@ -130,7 +217,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # Save a copy of the un-orthonormalized basis for later. # Convert it to ambient V (vector) coordinates while we're # at it, because we'd have to do it later anyway. - deortho_vector_basis = tuple( V(b.list()) for b in basis ) + deortho_vector_basis = tuple( V(_all2list(b)) for b in basis ) from mjo.eja.eja_utils import gram_schmidt basis = tuple(gram_schmidt(basis, inner_product)) @@ -142,7 +229,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # Now create the vector space for the algebra, which will have # its own set of non-ambient coordinates (in terms of the # supplied basis). - vector_basis = tuple( V(b.list()) for b in basis ) + vector_basis = tuple( V(_all2list(b)) for b in basis ) W = V.span_of_basis( vector_basis, check=check_axioms) if orthonormalize: @@ -174,7 +261,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # The jordan product returns a matrixy answer, so we # have to convert it to the algebra coordinates. elt = jordan_product(q_i, q_j) - elt = W.coordinate_vector(V(elt.list())) + elt = W.coordinate_vector(V(_all2list(elt))) self._multiplication_table[i][j] = self.from_vector(elt) if not orthonormalize: @@ -222,6 +309,35 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): def product_on_basis(self, i, j): + r""" + Returns the Jordan product of the `i` and `j`th basis elements. + + This completely defines the Jordan product on the algebra, and + is used direclty by our superclass machinery to implement + :meth:`product`. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS:: + + sage: set_random_seed() + sage: J = random_eja() + sage: n = J.dimension() + sage: ei = J.zero() + sage: ej = J.zero() + sage: ei_ej = J.zero()*J.zero() + sage: if n > 0: + ....: i = ZZ.random_element(n) + ....: j = ZZ.random_element(n) + ....: ei = J.gens()[i] + ....: ej = J.gens()[j] + ....: ei_ej = J.product_on_basis(i,j) + sage: ei*ej == ei_ej + True + + """ # We only stored the lower-triangular portion of the # multiplication table. if j <= i: @@ -279,22 +395,32 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: y = J.random_element() sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) True + """ B = self._inner_product_matrix return (B*x.to_vector()).inner_product(y.to_vector()) - def _is_commutative(self): + def is_associative(self): r""" - Whether or not this algebra's multiplication table is commutative. + Return whether or not this algebra's Jordan product is associative. + + SETUP:: + + sage: from mjo.eja.eja_algebra import ComplexHermitianEJA + + EXAMPLES:: + + sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False) + sage: J.is_associative() + False + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: A.is_associative() + True - 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.product_on_basis(i,j) == self.product_on_basis(i,j) - for i in range(self.dimension()) - for j in range(self.dimension()) ) + return "Associative" in self.category().axioms() def _is_jordanian(self): r""" @@ -303,13 +429,13 @@ 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.monomial(i)**2)*(self.monomial(i)*self.monomial(j)) + return all( (self.gens()[i]**2)*(self.gens()[i]*self.gens()[j]) == - (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j)) + (self.gens()[i])*((self.gens()[i]**2)*self.gens()[j]) for i in range(self.dimension()) for j in range(self.dimension()) ) @@ -322,26 +448,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.monomial(i) - y = self.monomial(j) - z = self.monomial(k) + x = self.gens()[i] + y = self.gens()[j] + z = self.gens()[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 @@ -377,6 +502,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): ... ValueError: not an element of this algebra + Tuples work as well, provided that the matrix basis for the + algebra consists of them:: + + sage: J1 = HadamardEJA(3) + 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) + TESTS: Ensure that we can convert any element of the two non-matrix @@ -393,13 +527,23 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: J(x.to_vector().column()) == x True + We cannot coerce elements between algebras just because their + matrix representations are compatible:: + + sage: J1 = HadamardEJA(3) + sage: J2 = JordanSpinEJA(3) + sage: J2(J1.one()) + Traceback (most recent call last): + ... + ValueError: not an element of this algebra + sage: J1(J2.zero()) + Traceback (most recent call last): + ... + ValueError: not an element of this algebra + """ msg = "not an element of this algebra" - if elt == 0: - # The superclass implementation of random_element() - # needs to be able to coerce "0" into the algebra. - return self.zero() - elif elt in self.base_ring(): + if elt in self.base_ring(): # Ensure that no base ring -> algebra coercion is performed # by this method. There's some stupidity in sage that would # otherwise propagate to this method; for example, sage thinks @@ -407,9 +551,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): raise ValueError(msg) try: + # Try to convert a vector into a column-matrix... elt = elt.column() except (AttributeError, TypeError): - # Try to convert a vector into a column-matrix + # and ignore failure, because we weren't really expecting + # a vector as an argument anyway. pass if elt not in self.matrix_space(): @@ -422,14 +568,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # closure whereas the base ring of the 3-by-3 identity matrix # could be QQ instead of QQbar. # + # And, we also have to handle Cartesian product bases (when + # the matrix basis consists of tuples) here. The "good news" + # is that we're already converting everything to long vectors, + # and that strategy works for tuples as well. + # # We pass check=False because the matrix basis is "guaranteed" # to be linearly independent... right? Ha ha. - V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols()) - W = V.span_of_basis( (_mat2vec(s) for s in self.matrix_basis()), + elt = _all2list(elt) + V = VectorSpace(self.base_ring(), len(elt)) + W = V.span_of_basis( (V(_all2list(s)) for s in self.matrix_basis()), check=False) try: - coords = W.coordinate_vector(_mat2vec(elt)) + coords = W.coordinate_vector(V(elt)) except ArithmeticError: # vector is not in free module raise ValueError(msg) @@ -656,8 +808,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # And to each subsequent row, prepend an entry that belongs to # the left-side "header column." - M += [ [self.monomial(i)] + [ self.product_on_basis(i,j) - for j in range(n) ] + M += [ [self.gens()[i]] + [ self.product_on_basis(i,j) + for j in range(n) ] for i in range(n) ] return table(M, header_row=True, header_column=True, frame=True) @@ -687,7 +839,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Why implement this for non-matrix algebras? Avoiding special cases for the :class:`BilinearFormEJA` pays with simplicity in its own right. But mainly, we would like to be able to assume - that elements of a :class:`DirectSumEJA` can be displayed + that elements of a :class:`CartesianProductEJA` can be displayed nicely, without having to have special classes for direct sums one of whose components was a matrix algebra. @@ -727,12 +879,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) @@ -993,14 +1182,12 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): if not c.is_idempotent(): raise ValueError("element is not idempotent: %s" % c) - from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra - # Default these to what they should be if they turn out to be # trivial, because eigenspaces_left() won't return eigenvalues # corresponding to trivial spaces (e.g. it returns only the # eigenspace corresponding to lambda=1 if you take the # decomposition relative to the identity element). - trivial = FiniteDimensionalEJASubalgebra(self, ()) + trivial = self.subalgebra(()) J0 = trivial # eigenvalue zero J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half J1 = trivial # eigenvalue one @@ -1010,9 +1197,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): J5 = eigspace else: gens = tuple( self.from_vector(b) for b in eigspace.basis() ) - subalg = FiniteDimensionalEJASubalgebra(self, - gens, - check_axioms=False) + subalg = self.subalgebra(gens, check_axioms=False) if eigval == 0: J0 = subalg elif eigval == 1: @@ -1125,7 +1310,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.monomial(k).operator().matrix()[i,j] + return sum( vars[k]*self.gens()[k].operator().matrix()[i,j] for k in range(n) ) L_x = matrix(F, n, n, L_x_i_j) @@ -1151,10 +1336,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # The theory says that only the first "r" coefficients are # nonzero, and they actually live in the original polynomial - # ring and not the fraction field. We negate them because - # in the actual characteristic polynomial, they get moved - # to the other side where x^r lives. - return -A_rref.solve_right(E*b).change_ring(R)[:r] + # ring and not the fraction field. We negate them because in + # the actual characteristic polynomial, they get moved to the + # other side where x^r lives. We don't bother to trim A_rref + # down to a square matrix and solve the resulting system, + # because the upper-left r-by-r portion of A_rref is + # guaranteed to be the identity matrix, so e.g. + # + # A_rref.solve_right(Y) + # + # would just be returning Y. + return (-E*b)[:r].change_ring(R) @cached_method def rank(self): @@ -1215,7 +1407,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: set_random_seed() # long time sage: J = random_eja() # long time - sage: caches = J.rank() # long time + sage: cached = J.rank() # long time sage: J.rank.clear_cache() # long time sage: J.rank() == cached # long time True @@ -1224,6 +1416,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): return len(self._charpoly_coefficients()) + def subalgebra(self, basis, **kwargs): + r""" + Create a subalgebra of this algebra from the given basis. + """ + from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra + return FiniteDimensionalEJASubalgebra(self, basis, **kwargs) + + def vector_space(self): """ Return the vector space that underlies this algebra. @@ -1242,7 +1442,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): return self.zero().to_vector().parent().ambient_vector_space() - Element = FiniteDimensionalEJAElement class RationalBasisEJA(FiniteDimensionalEJA): r""" @@ -1576,9 +1775,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 @@ -1844,9 +2043,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 @@ -2141,9 +2340,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 @@ -2351,7 +2550,11 @@ class HadamardEJA(ConcreteEJA): if "check_axioms" not in kwargs: kwargs["check_axioms"] = False column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() ) - super().__init__(column_basis, jordan_product, inner_product, **kwargs) + super().__init__(column_basis, + jordan_product, + inner_product, + associative=True, + **kwargs) self.rank.set_cache(n) if n == 0: @@ -2668,262 +2871,432 @@ class TrivialEJA(ConcreteEJA): return cls(**kwargs) -class DirectSumEJA(FiniteDimensionalEJA): +class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct, + FiniteDimensionalEJA): r""" - The external (orthogonal) direct sum of two other Euclidean Jordan - algebras. Essentially the Cartesian product of its two factors. - Every Euclidean Jordan algebra decomposes into an orthogonal - direct sum of simple Euclidean Jordan algebras, so no generality - is lost by providing only this construction. + The external (orthogonal) direct sum of two or more Euclidean + Jordan algebras. Every Euclidean Jordan algebra decomposes into an + orthogonal direct sum of simple Euclidean Jordan algebras which is + then isometric to a Cartesian product, so no generality is lost by + providing only this construction. SETUP:: sage: from mjo.eja.eja_algebra import (random_eja, + ....: CartesianProductEJA, ....: HadamardEJA, - ....: RealSymmetricEJA, - ....: DirectSumEJA) + ....: JordanSpinEJA, + ....: RealSymmetricEJA) - EXAMPLES:: + EXAMPLES: + + The Jordan product is inherited from our factors and implemented by + our CombinatorialFreeModule Cartesian product superclass:: + sage: set_random_seed() sage: J1 = HadamardEJA(2) - sage: J2 = RealSymmetricEJA(3) - sage: J = DirectSumEJA(J1,J2) - sage: J.dimension() - 8 + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: x,y = J.random_elements(2) + sage: x*y in J + True + + The ability to retrieve the original factors is implemented by our + CombinatorialFreeModule Cartesian product superclass:: + + sage: J1 = HadamardEJA(2, field=QQ) + sage: J2 = JordanSpinEJA(3, field=QQ) + sage: J = cartesian_product([J1,J2]) + sage: J.cartesian_factors() + (Euclidean Jordan algebra of dimension 2 over Rational Field, + Euclidean Jordan algebra of dimension 3 over Rational Field) + + You can provide more than two factors:: + + sage: J1 = HadamardEJA(2) + sage: J2 = JordanSpinEJA(3) + sage: J3 = RealSymmetricEJA(3) + sage: cartesian_product([J1,J2,J3]) + Euclidean Jordan algebra of dimension 2 over Algebraic Real + Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic + Real Field (+) Euclidean Jordan algebra of dimension 6 over + Algebraic Real Field + + Rank is additive on a Cartesian product:: + + sage: J1 = HadamardEJA(1) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J1.rank.clear_cache() + sage: J2.rank.clear_cache() + sage: J.rank.clear_cache() sage: J.rank() - 5 - sage: J.matrix_space() - The Cartesian product of (Full MatrixSpace of 2 by 1 dense matrices - over Algebraic Real Field, Full MatrixSpace of 3 by 3 dense matrices - over Algebraic Real Field) + 3 + sage: J.rank() == J1.rank() + J2.rank() + True + + The same rank computation works over the rationals, with whatever + basis you like:: + + sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False) + sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False) + sage: J = cartesian_product([J1,J2]) + sage: J1.rank.clear_cache() + sage: J2.rank.clear_cache() + sage: J.rank.clear_cache() + sage: J.rank() + 3 + sage: J.rank() == J1.rank() + J2.rank() + True + + The product algebra will be associative if and only if all of its + components are associative:: + + sage: J1 = HadamardEJA(2) + sage: J1.is_associative() + True + sage: J2 = HadamardEJA(3) + sage: J2.is_associative() + True + sage: J3 = RealSymmetricEJA(3) + sage: J3.is_associative() + False + sage: CP1 = cartesian_product([J1,J2]) + sage: CP1.is_associative() + True + sage: CP2 = cartesian_product([J1,J3]) + sage: CP2.is_associative() + False TESTS: - The external direct sum construction is only valid when the two factors - have the same base ring; an error is raised otherwise:: + All factors must share the same base field:: - sage: set_random_seed() - sage: J1 = random_eja(field=AA) - sage: J2 = random_eja(field=QQ,orthonormalize=False) - sage: J = DirectSumEJA(J1,J2) + sage: J1 = HadamardEJA(2, field=QQ) + sage: J2 = RealSymmetricEJA(2) + sage: CartesianProductEJA((J1,J2)) Traceback (most recent call last): ... - ValueError: algebras must share the same base field + ValueError: all factors must share the same base field + + The cached unit element is the same one that would be computed:: + + sage: set_random_seed() # long time + sage: J1 = random_eja() # long time + sage: J2 = random_eja() # long time + sage: J = cartesian_product([J1,J2]) # long time + sage: actual = J.one() # long time + sage: J.one.clear_cache() # long time + sage: expected = J.one() # long time + sage: actual == expected # long time + True """ - def __init__(self, J1, J2, **kwargs): - if J1.base_ring() != J2.base_ring(): - raise ValueError("algebras must share the same base field") - field = J1.base_ring() - - M = J1.matrix_space().cartesian_product(J2.matrix_space()) - self._cartprod_algebra = J1.cartesian_product(J2) - - self._matrix_basis = tuple( [M((a,0)) for a in J1.matrix_basis()] + - [M((0,b)) for b in J2.matrix_basis()] ) - - n = len(self._matrix_basis) - self._sets = None - CombinatorialFreeModule.__init__( - self, - field, - range(n), - category=self._cartprod_algebra.category(), - bracket=False, - **kwargs) - self.rank.set_cache(J1.rank() + J2.rank()) + 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 ): + raise ValueError("all factors must share the same base field") + + associative = all( m.is_associative() for m in algebras ) + + # 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() + ) - def product(self,x,y): + # 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) + )) + + def inner_product(x, y): + return sum( + Js[i](x[i]).inner_product(Js[i](y[i])) for i in range(m) + ) + + # There's no need to check the field since it already came + # from an EJA. Likewise the axioms are guaranteed to be + # satisfied, unless the guy writing this class sucks. + # + # If you want the basis to be orthonormalized, orthonormalize + # the factors. + FiniteDimensionalEJA.__init__(self, + basis, + jordan_product, + inner_product, + field=field, + orthonormalize=False, + associative=associative, + cartesian_product=True, + 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)) + + def matrix_space(self): r""" + Return the space that our matrix basis lives in as a Cartesian + product. + SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, - ....: ComplexHermitianEJA, - ....: DirectSumEJA) + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: RealSymmetricEJA) - TESTS:: + EXAMPLES:: + + sage: J1 = HadamardEJA(1) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J.matrix_space() + The Cartesian product of (Full MatrixSpace of 1 by 1 dense + matrices over Algebraic Real Field, Full MatrixSpace of 2 + by 2 dense matrices over Algebraic Real Field) + + """ + from sage.categories.cartesian_product import cartesian_product + return cartesian_product( [J.matrix_space() + for J in self.cartesian_factors()] ) + + @cached_method + def cartesian_projection(self, i): + r""" + SETUP:: + + sage: from mjo.eja.eja_algebra import (random_eja, + ....: JordanSpinEJA, + ....: HadamardEJA, + ....: RealSymmetricEJA, + ....: ComplexHermitianEJA) + + EXAMPLES: + + The projection morphisms are Euclidean Jordan algebra + operators:: + + sage: J1 = HadamardEJA(2) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J.cartesian_projection(0) + Linear operator between finite-dimensional Euclidean Jordan + algebras represented by the matrix: + [1 0 0 0 0] + [0 1 0 0 0] + Domain: Euclidean Jordan algebra of dimension 2 over Algebraic + Real Field (+) Euclidean Jordan algebra of dimension 3 over + Algebraic Real Field + Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic + Real Field + sage: J.cartesian_projection(1) + Linear operator between finite-dimensional Euclidean Jordan + algebras represented by the matrix: + [0 0 1 0 0] + [0 0 0 1 0] + [0 0 0 0 1] + Domain: Euclidean Jordan algebra of dimension 2 over Algebraic + Real Field (+) Euclidean Jordan algebra of dimension 3 over + Algebraic Real Field + Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic + Real Field + + The projections work the way you'd expect on the vector + representation of an element:: + + sage: J1 = JordanSpinEJA(2) + sage: J2 = ComplexHermitianEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: pi_left = J.cartesian_projection(0) + sage: pi_right = J.cartesian_projection(1) + sage: pi_left(J.one()).to_vector() + (1, 0) + sage: pi_right(J.one()).to_vector() + (1, 0, 0, 1) + sage: J.one().to_vector() + (1, 0, 1, 0, 0, 1) + + TESTS: + + The answer never changes:: sage: set_random_seed() - sage: J1 = JordanSpinEJA(3, field=QQ) - sage: J2 = ComplexHermitianEJA(2, field=QQ, orthonormalize=False) - sage: J = DirectSumEJA(J1,J2) - sage: J.random_element()*J.random_element() in J + sage: J1 = random_eja() + sage: J2 = random_eja() + sage: J = cartesian_product([J1,J2]) + sage: P0 = J.cartesian_projection(0) + sage: P1 = J.cartesian_projection(0) + sage: P0 == P1 True """ - xv = self._cartprod_algebra.from_vector(x.to_vector()) - yv = self._cartprod_algebra.from_vector(y.to_vector()) - return self.from_vector((xv*yv).to_vector()) - + Ji = self.cartesian_factors()[i] + # Requires the fix on Trac 31421/31422 to work! + Pi = super().cartesian_projection(i) + return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix()) - def cartesian_factors(self): + @cached_method + def cartesian_embedding(self, i): r""" - Return the pair of this algebra's factors. - SETUP:: - sage: from mjo.eja.eja_algebra import (HadamardEJA, + sage: from mjo.eja.eja_algebra import (random_eja, ....: JordanSpinEJA, - ....: DirectSumEJA) + ....: HadamardEJA, + ....: RealSymmetricEJA) - EXAMPLES:: + EXAMPLES: + + The embedding morphisms are Euclidean Jordan algebra + operators:: + + sage: J1 = HadamardEJA(2) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J.cartesian_embedding(0) + Linear operator between finite-dimensional Euclidean Jordan + algebras represented by the matrix: + [1 0] + [0 1] + [0 0] + [0 0] + [0 0] + Domain: Euclidean Jordan algebra of dimension 2 over + Algebraic Real Field + Codomain: Euclidean Jordan algebra of dimension 2 over + Algebraic Real Field (+) Euclidean Jordan algebra of + dimension 3 over Algebraic Real Field + sage: J.cartesian_embedding(1) + Linear operator between finite-dimensional Euclidean Jordan + algebras represented by the matrix: + [0 0 0] + [0 0 0] + [1 0 0] + [0 1 0] + [0 0 1] + Domain: Euclidean Jordan algebra of dimension 3 over + Algebraic Real Field + Codomain: Euclidean Jordan algebra of dimension 2 over + Algebraic Real Field (+) Euclidean Jordan algebra of + dimension 3 over Algebraic Real Field + + The embeddings work the way you'd expect on the vector + representation of an element:: + + sage: J1 = JordanSpinEJA(3) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: iota_left = J.cartesian_embedding(0) + sage: iota_right = J.cartesian_embedding(1) + sage: iota_left(J1.zero()) == J.zero() + True + sage: iota_right(J2.zero()) == J.zero() + True + sage: J1.one().to_vector() + (1, 0, 0) + sage: iota_left(J1.one()).to_vector() + (1, 0, 0, 0, 0, 0) + sage: J2.one().to_vector() + (1, 0, 1) + sage: iota_right(J2.one()).to_vector() + (0, 0, 0, 1, 0, 1) + sage: J.one().to_vector() + (1, 0, 0, 1, 0, 1) + + TESTS: + + The answer never changes:: + + sage: set_random_seed() + sage: J1 = random_eja() + sage: J2 = random_eja() + sage: J = cartesian_product([J1,J2]) + sage: E0 = J.cartesian_embedding(0) + sage: E1 = J.cartesian_embedding(0) + sage: E0 == E1 + True + + Composing a projection with the corresponding inclusion should + produce the identity map, and mismatching them should produce + the zero map:: + + sage: set_random_seed() + sage: J1 = random_eja() + sage: J2 = random_eja() + sage: J = cartesian_product([J1,J2]) + sage: iota_left = J.cartesian_embedding(0) + sage: iota_right = J.cartesian_embedding(1) + sage: pi_left = J.cartesian_projection(0) + sage: pi_right = J.cartesian_projection(1) + sage: pi_left*iota_left == J1.one().operator() + True + sage: pi_right*iota_right == J2.one().operator() + True + sage: (pi_left*iota_right).is_zero() + True + sage: (pi_right*iota_left).is_zero() + True + + """ + Ji = self.cartesian_factors()[i] + # Requires the fix on Trac 31421/31422 to work! + Ei = super().cartesian_embedding(i) + return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix()) + + + +FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA + +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) - sage: J1 = HadamardEJA(2, field=QQ) - sage: J2 = JordanSpinEJA(3, field=QQ) - sage: J = DirectSumEJA(J1,J2) - sage: J.cartesian_factors() - (Euclidean Jordan algebra of dimension 2 over Rational Field, - Euclidean Jordan algebra of dimension 3 over Rational Field) - - """ - return self._cartprod_algebra.cartesian_factors() - - -# def projections(self): -# r""" -# Return a pair of projections onto this algebra's factors. - -# SETUP:: - -# sage: from mjo.eja.eja_algebra import (JordanSpinEJA, -# ....: ComplexHermitianEJA, -# ....: DirectSumEJA) - -# EXAMPLES:: - -# sage: J1 = JordanSpinEJA(2) -# sage: J2 = ComplexHermitianEJA(2) -# sage: J = DirectSumEJA(J1,J2) -# sage: (pi_left, pi_right) = J.projections() -# sage: J.one().to_vector() -# (1, 0, 1, 0, 0, 1) -# sage: pi_left(J.one()).to_vector() -# (1, 0) -# sage: pi_right(J.one()).to_vector() -# (1, 0, 0, 1) - -# """ -# (J1,J2) = self.factors() -# m = J1.dimension() -# n = J2.dimension() -# V_basis = self.vector_space().basis() -# # Need to specify the dimensions explicitly so that we don't -# # wind up with a zero-by-zero matrix when we want e.g. a -# # zero-by-two matrix (important for composing things). -# P1 = matrix(self.base_ring(), m, m+n, V_basis[:m]) -# P2 = matrix(self.base_ring(), n, m+n, V_basis[m:]) -# pi_left = FiniteDimensionalEJAOperator(self,J1,P1) -# pi_right = FiniteDimensionalEJAOperator(self,J2,P2) -# return (pi_left, pi_right) - -# def inclusions(self): -# r""" -# Return the pair of inclusion maps from our factors into us. - -# SETUP:: - -# sage: from mjo.eja.eja_algebra import (random_eja, -# ....: JordanSpinEJA, -# ....: RealSymmetricEJA, -# ....: DirectSumEJA) - -# EXAMPLES:: - -# sage: J1 = JordanSpinEJA(3) -# sage: J2 = RealSymmetricEJA(2) -# sage: J = DirectSumEJA(J1,J2) -# sage: (iota_left, iota_right) = J.inclusions() -# sage: iota_left(J1.zero()) == J.zero() -# True -# sage: iota_right(J2.zero()) == J.zero() -# True -# sage: J1.one().to_vector() -# (1, 0, 0) -# sage: iota_left(J1.one()).to_vector() -# (1, 0, 0, 0, 0, 0) -# sage: J2.one().to_vector() -# (1, 0, 1) -# sage: iota_right(J2.one()).to_vector() -# (0, 0, 0, 1, 0, 1) -# sage: J.one().to_vector() -# (1, 0, 0, 1, 0, 1) - -# TESTS: - -# Composing a projection with the corresponding inclusion should -# produce the identity map, and mismatching them should produce -# the zero map:: - -# sage: set_random_seed() -# sage: J1 = random_eja() -# sage: J2 = random_eja() -# sage: J = DirectSumEJA(J1,J2) -# sage: (iota_left, iota_right) = J.inclusions() -# sage: (pi_left, pi_right) = J.projections() -# sage: pi_left*iota_left == J1.one().operator() -# True -# sage: pi_right*iota_right == J2.one().operator() -# True -# sage: (pi_left*iota_right).is_zero() -# True -# sage: (pi_right*iota_left).is_zero() -# True - -# """ -# (J1,J2) = self.factors() -# m = J1.dimension() -# n = J2.dimension() -# V_basis = self.vector_space().basis() -# # Need to specify the dimensions explicitly so that we don't -# # wind up with a zero-by-zero matrix when we want e.g. a -# # two-by-zero matrix (important for composing things). -# I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m]) -# I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:]) -# iota_left = FiniteDimensionalEJAOperator(J1,self,I1) -# iota_right = FiniteDimensionalEJAOperator(J2,self,I2) -# return (iota_left, iota_right) - -# def inner_product(self, x, y): -# r""" -# The standard Cartesian inner-product. - -# We project ``x`` and ``y`` onto our factors, and add up the -# inner-products from the subalgebras. - -# SETUP:: - - -# sage: from mjo.eja.eja_algebra import (HadamardEJA, -# ....: QuaternionHermitianEJA, -# ....: DirectSumEJA) - -# EXAMPLE:: - -# sage: J1 = HadamardEJA(3,field=QQ) -# sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False) -# sage: J = DirectSumEJA(J1,J2) -# sage: x1 = J1.one() -# sage: x2 = x1 -# sage: y1 = J2.one() -# sage: y2 = y1 -# sage: x1.inner_product(x2) -# 3 -# sage: y1.inner_product(y2) -# 2 -# sage: J.one().inner_product(J.one()) -# 5 - -# """ -# (pi_left, pi_right) = self.projections() -# x1 = pi_left(x) -# x2 = pi_right(x) -# y1 = pi_left(y) -# y2 = pi_right(y) - -# return (x1.inner_product(y1) + x2.inner_product(y2)) + 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 random_eja = ConcreteEJA.random_instance