X-Git-Url: http://gitweb.michael.orlitzky.com/?p=sage.d.git;a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=5bf597565b2ae292032c3f0a532763d7889cd2a8;hp=db899dd984594bf8ef88cb0e44f3812d2e843996;hb=db1f7761ebf564221669137ae07476ea45d82a2c;hpb=3b3755e21ddd4f1c44fd8a6bb179f4c3a59c5d2d diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index db899dd..5bf5975 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,13 +57,13 @@ 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.categories.sets_cat import cartesian_product from sage.combinat.free_module import CombinatorialFreeModule from sage.matrix.constructor import matrix from sage.matrix.matrix_space import MatrixSpace @@ -31,7 +75,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 +83,50 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): INPUT: - - basis -- a tuple of basis elements in their matrix form. + - ``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: - - 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. + We should compute that an element subalgebra is associative even + if we circumvent the element method:: - - 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. + 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 @@ -59,10 +137,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'): + prefix="b"): + + n = len(basis) if check_field: if not field.is_subring(RR): @@ -71,10 +152,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # we've specified a real embedding. raise ValueError("scalar field is not real") - # 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 check_axioms: # Check commutativity of the Jordan and inner-products. # This has to be done before we build the multiplication @@ -92,31 +169,48 @@ 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: + # 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. - 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, # as well as a subspace W of V spanned by those (vectorized) # basis elements. The W-coordinates are the coefficients that - # we see in things like x = 1*e1 + 2*e2. + # we see in things like x = 1*b1 + 2*b2. vector_basis = basis 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 +224,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 +236,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 +268,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 +316,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: bi = J.zero() + sage: bj = J.zero() + sage: bi_bj = J.zero()*J.zero() + sage: if n > 0: + ....: i = ZZ.random_element(n) + ....: j = ZZ.random_element(n) + ....: bi = J.monomial(i) + ....: bj = J.monomial(j) + ....: bi_bj = J.product_on_basis(i,j) + sage: bi*bj == bi_bj + True + + """ # We only stored the lower-triangular portion of the # multiplication table. if j <= i: @@ -279,11 +402,33 @@ 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_associative(self): + r""" + 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 + + """ + return "Associative" in self.category().axioms() + def _is_commutative(self): r""" Whether or not this algebra's multiplication table is commutative. @@ -292,9 +437,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): 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 all( x*y == y*x for x in self.gens() for y in self.gens() ) def _is_jordanian(self): r""" @@ -313,6 +456,92 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): 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 @@ -322,11 +551,14 @@ 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()): @@ -336,12 +568,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): 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 @@ -355,7 +583,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) @@ -377,29 +606,42 @@ 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]) ) + b1 + b5 + 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 + 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 +649,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 +666,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) @@ -638,15 +888,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: J = JordanSpinEJA(4) sage: J.multiplication_table() +----++----+----+----+----+ - | * || e0 | e1 | e2 | e3 | + | * || b0 | b1 | b2 | b3 | +====++====+====+====+====+ - | e0 || e0 | e1 | e2 | e3 | + | b0 || b0 | b1 | b2 | b3 | +----++----+----+----+----+ - | e1 || e1 | e0 | 0 | 0 | + | b1 || b1 | b0 | 0 | 0 | +----++----+----+----+----+ - | e2 || e2 | 0 | e0 | 0 | + | b2 || b2 | 0 | b0 | 0 | +----++----+----+----+----+ - | e3 || e3 | 0 | 0 | e0 | + | b3 || b3 | 0 | 0 | b0 | +----++----+----+----+----+ """ @@ -656,8 +906,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.monomial(i)] + [ self.monomial(i)*self.monomial(j) + for j in range(n) ] for i in range(n) ] return table(M, header_row=True, header_column=True, frame=True) @@ -687,7 +937,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. @@ -700,7 +950,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: J = RealSymmetricEJA(2) sage: J.basis() - Finite family {0: e0, 1: e1, 2: e2} + Finite family {0: b0, 1: b1, 2: b2} sage: J.matrix_basis() ( [1 0] [ 0 0.7071067811865475?] [0 0] @@ -711,7 +961,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: J = JordanSpinEJA(2) sage: J.basis() - Finite family {0: e0, 1: e1} + Finite family {0: b0, 1: b1} sage: J.matrix_basis() ( [1] [0] @@ -727,17 +977,54 @@ 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) else: - return self._matrix_basis[0].matrix_space() + return self.matrix_basis()[0].parent() @cached_method @@ -756,20 +1043,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: J = HadamardEJA(5) sage: J.one() - e0 + e1 + e2 + e3 + e4 + b0 + b1 + b2 + b3 + b4 The unit element in the Hadamard EJA is inherited in the subalgebras generated by its elements:: sage: J = HadamardEJA(5) sage: J.one() - e0 + e1 + e2 + e3 + e4 + b0 + b1 + b2 + b3 + b4 sage: x = sum(J.gens()) sage: A = x.subalgebra_generated_by(orthonormalize=False) sage: A.one() - f0 + c0 sage: A.one().superalgebra_element() - e0 + e1 + e2 + e3 + e4 + b0 + b1 + b2 + b3 + b4 TESTS: @@ -993,14 +1280,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 +1295,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: @@ -1101,6 +1384,21 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): r""" The `r` polynomial coefficients of the "characteristic polynomial of" function. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS: + + The theory shows that these are all homogeneous polynomials of + a known degree:: + + sage: set_random_seed() + sage: J = random_eja() + sage: all(p.is_homogeneous() for p in J._charpoly_coefficients()) + True + """ n = self.dimension() R = self.coordinate_polynomial_ring() @@ -1136,10 +1434,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): @@ -1200,7 +1505,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 @@ -1209,6 +1514,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. @@ -1227,7 +1540,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): return self.zero().to_vector().parent().ambient_vector_space() - Element = FiniteDimensionalEJAElement class RationalBasisEJA(FiniteDimensionalEJA): r""" @@ -1265,6 +1577,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 @@ -1278,17 +1597,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""" @@ -1412,6 +1725,21 @@ class ConcreteEJA(RationalBasisEJA): class MatrixEJA: + @staticmethod + def jordan_product(X,Y): + return (X*Y + Y*X)/2 + + @staticmethod + def trace_inner_product(X,Y): + r""" + A trace inner-product for matrices that aren't embedded in the + reals. + """ + # We take the norm (absolute value) because Octonions() isn't + # smart enough yet to coerce its one() into the base field. + return (X*Y).trace().abs() + +class RealEmbeddedMatrixEJA(MatrixEJA): @staticmethod def dimension_over_reals(): r""" @@ -1457,9 +1785,6 @@ class MatrixEJA: raise ValueError("the matrix 'M' must be a real embedding") return M - @staticmethod - def jordan_product(X,Y): - return (X*Y + Y*X)/2 @classmethod def trace_inner_product(cls,X,Y): @@ -1468,29 +1793,11 @@ class MatrixEJA: SETUP:: - sage: from mjo.eja.eja_algebra import (RealSymmetricEJA, - ....: ComplexHermitianEJA, + sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, ....: QuaternionHermitianEJA) EXAMPLES:: - This gives the same answer as it would if we computed the trace - from the unembedded (original) matrices:: - - sage: set_random_seed() - sage: J = RealSymmetricEJA.random_instance() - sage: x,y = J.random_elements(2) - sage: Xe = x.to_matrix() - sage: Ye = y.to_matrix() - sage: X = J.real_unembed(Xe) - sage: Y = J.real_unembed(Ye) - sage: expected = (X*Y).trace() - sage: actual = J.trace_inner_product(Xe,Ye) - sage: actual == expected - True - - :: - sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() sage: x,y = J.random_elements(2) @@ -1518,27 +1825,15 @@ class MatrixEJA: True """ - Xu = cls.real_unembed(X) - Yu = cls.real_unembed(Y) - tr = (Xu*Yu).trace() + # This does in fact compute the real part of the trace. + # If we compute the trace of e.g. a complex matrix M, + # then we do so by adding up its diagonal entries -- + # call them z_1 through z_n. The real embedding of z_1 + # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace + # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth. + return (X*Y).trace()/cls.dimension_over_reals() - try: - # Works in QQ, AA, RDF, et cetera. - return tr.real() - except AttributeError: - # A quaternion doesn't have a real() method, but does - # have coefficient_tuple() method that returns the - # coefficients of 1, i, j, and k -- in that order. - return tr.coefficient_tuple()[0] - - -class RealMatrixEJA(MatrixEJA): - @staticmethod - def dimension_over_reals(): - return 1 - - -class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): +class RealSymmetricEJA(ConcreteEJA, MatrixEJA): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -1551,19 +1846,19 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): EXAMPLES:: sage: J = RealSymmetricEJA(2) - sage: e0, e1, e2 = J.gens() - sage: e0*e0 - e0 - sage: e1*e1 - 1/2*e0 + 1/2*e2 - sage: e2*e2 - e2 + sage: b0, b1, b2 = J.gens() + sage: b0*b0 + b0 + sage: b1*b1 + 1/2*b0 + 1/2*b2 + sage: b2*b2 + b2 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 @@ -1604,7 +1899,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): """ @classmethod - def _denormalized_basis(cls, n): + def _denormalized_basis(cls, n, field): """ Return a basis for the space of real symmetric n-by-n matrices. @@ -1616,7 +1911,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): sage: set_random_seed() sage: n = ZZ.random_element(1,5) - sage: B = RealSymmetricEJA._denormalized_basis(n) + sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ) sage: all( M.is_symmetric() for M in B) True @@ -1626,7 +1921,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): S = [] for i in range(n): for j in range(i+1): - Eij = matrix(ZZ, n, lambda k,l: k==i and l==j) + Eij = matrix(field, n, lambda k,l: k==i and l==j) if i == j: Sij = Eij else: @@ -1647,26 +1942,32 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, **kwargs) - def __init__(self, n, **kwargs): + def __init__(self, n, field=AA, **kwargs): # We know this is a valid EJA, but will double-check # if the user passes check_axioms=True. 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,field), + self.jordan_product, + self.trace_inner_product, + field=field, + 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(). self.rank.set_cache(n) - idV = matrix.identity(ZZ, self.dimension_over_reals()*n) + idV = self.matrix_space().one() self.one.set_cache(self(idV)) -class ComplexMatrixEJA(MatrixEJA): +class ComplexMatrixEJA(RealEmbeddedMatrixEJA): # A manual dictionary-cache for the complex_extension() method, # since apparently @classmethods can't also be @cached_methods. _complex_extension = {} @@ -1745,7 +2046,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... @@ -1792,7 +2093,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()) @@ -1829,9 +2130,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 @@ -1873,7 +2174,7 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): """ @classmethod - def _denormalized_basis(cls, n): + def _denormalized_basis(cls, n, field): """ Returns a basis for the space of complex Hermitian n-by-n matrices. @@ -1891,15 +2192,14 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): sage: set_random_seed() sage: n = ZZ.random_element(1,5) - sage: B = ComplexHermitianEJA._denormalized_basis(n) + sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ) sage: all( M.is_symmetric() for M in B) True """ - field = ZZ - R = PolynomialRing(field, 'z') + R = PolynomialRing(ZZ, 'z') z = R.gen() - F = field.extension(z**2 + 1, 'I') + F = ZZ.extension(z**2 + 1, 'I') I = F.gen(1) # This is like the symmetric case, but we need to be careful: @@ -1930,20 +2230,26 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): # "erase" E_ij Eij[i,j] = 0 - # Since we embedded these, we can drop back to the "field" that we - # started with instead of the complex extension "F". + # Since we embedded the entries, we can drop back to the + # desired real "field" instead of the extension "F". return tuple( s.change_ring(field) for s in S ) - def __init__(self, n, **kwargs): + def __init__(self, n, field=AA, **kwargs): # We know this is a valid EJA, but will double-check # if the user passes check_axioms=True. 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,field), + self.jordan_product, + self.trace_inner_product, + field=field, + 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(). @@ -1963,7 +2269,7 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, **kwargs) -class QuaternionMatrixEJA(MatrixEJA): +class QuaternionMatrixEJA(RealEmbeddedMatrixEJA): # A manual dictionary-cache for the quaternion_extension() method, # since apparently @classmethods can't also be @cached_methods. @@ -2026,7 +2332,7 @@ class QuaternionMatrixEJA(MatrixEJA): True """ - super(QuaternionMatrixEJA,cls).real_embed(M) + super().real_embed(M) quaternions = M.base_ring() n = M.nrows() @@ -2081,7 +2387,7 @@ class QuaternionMatrixEJA(MatrixEJA): True """ - super(QuaternionMatrixEJA,cls).real_unembed(M) + super().real_unembed(M) n = ZZ(M.nrows()) d = cls.dimension_over_reals() @@ -2126,9 +2432,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 @@ -2169,7 +2475,7 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA): """ @classmethod - def _denormalized_basis(cls, n): + def _denormalized_basis(cls, n, field): """ Returns a basis for the space of quaternion Hermitian n-by-n matrices. @@ -2187,12 +2493,11 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA): sage: set_random_seed() sage: n = ZZ.random_element(1,5) - sage: B = QuaternionHermitianEJA._denormalized_basis(n) + sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ) sage: all( M.is_symmetric() for M in B ) True """ - field = ZZ Q = QuaternionAlgebra(QQ,-1,-1) I,J,K = Q.gens() @@ -2236,20 +2541,27 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA): # "erase" E_ij Eij[i,j] = 0 - # Since we embedded these, we can drop back to the "field" that we - # started with instead of the quaternion algebra "Q". + # Since we embedded the entries, we can drop back to the + # desired real "field" instead of the quaternion algebra "Q". return tuple( s.change_ring(field) for s in S ) - def __init__(self, n, **kwargs): + def __init__(self, n, field=AA, **kwargs): # We know this is a valid EJA, but will double-check # if the user passes check_axioms=True. 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,field), + self.jordan_product, + self.trace_inner_product, + field=field, + 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(). @@ -2292,19 +2604,19 @@ class HadamardEJA(ConcreteEJA): This multiplication table can be verified by hand:: sage: J = HadamardEJA(3) - sage: e0,e1,e2 = J.gens() - sage: e0*e0 - e0 - sage: e0*e1 + sage: b0,b1,b2 = J.gens() + sage: b0*b0 + b0 + sage: b0*b1 0 - sage: e0*e2 + sage: b0*b2 0 - sage: e1*e1 - e1 - sage: e1*e2 + sage: b1*b1 + b1 + sage: b1*b2 0 - sage: e2*e2 - e2 + sage: b2*b2 + b2 TESTS: @@ -2314,7 +2626,7 @@ class HadamardEJA(ConcreteEJA): (r0, r1, r2) """ - def __init__(self, n, **kwargs): + def __init__(self, n, field=AA, **kwargs): if n == 0: jordan_product = lambda x,y: x inner_product = lambda x,y: x @@ -2335,8 +2647,14 @@ class HadamardEJA(ConcreteEJA): if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False 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) + column_basis = tuple( b.column() + for b in FreeModule(field, n).basis() ) + super().__init__(column_basis, + jordan_product, + inner_product, + field=field, + associative=True, + **kwargs) self.rank.set_cache(n) if n == 0: @@ -2442,7 +2760,7 @@ class BilinearFormEJA(ConcreteEJA): True """ - def __init__(self, B, **kwargs): + def __init__(self, B, field=AA, **kwargs): # The matrix "B" is supplied by the user in most cases, # so it makes sense to check whether or not its positive- # definite unless we are specifically asked not to... @@ -2470,11 +2788,20 @@ class BilinearFormEJA(ConcreteEJA): return P([z0] + zbar.list()) 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) + column_basis = tuple( b.column() + for b in FreeModule(field, n).basis() ) + + # 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, + field=field, + associative=associative, + **kwargs) # The rank of this algebra is two, unless we're in a # one-dimensional ambient space (because the rank is bounded @@ -2534,20 +2861,20 @@ class JordanSpinEJA(BilinearFormEJA): This multiplication table can be verified by hand:: sage: J = JordanSpinEJA(4) - sage: e0,e1,e2,e3 = J.gens() - sage: e0*e0 - e0 - sage: e0*e1 - e1 - sage: e0*e2 - e2 - sage: e0*e3 - e3 - sage: e1*e2 + sage: b0,b1,b2,b3 = J.gens() + sage: b0*b0 + b0 + sage: b0*b1 + b1 + sage: b0*b2 + b2 + sage: b0*b3 + b3 + sage: b1*b2 0 - sage: e1*e3 + sage: b1*b3 0 - sage: e2*e3 + sage: b2*b3 0 We can change the generator prefix:: @@ -2568,7 +2895,7 @@ class JordanSpinEJA(BilinearFormEJA): True """ - def __init__(self, n, **kwargs): + def __init__(self, n, *args, **kwargs): # This is a special case of the BilinearFormEJA with the # identity matrix as its bilinear form. B = matrix.identity(ZZ, n) @@ -2579,7 +2906,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, *args, **kwargs) @staticmethod def _max_random_instance_size(): @@ -2637,10 +2964,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) @@ -2652,228 +2981,495 @@ class TrivialEJA(ConcreteEJA): # inappropriate for us. return cls(**kwargs) -class DirectSumEJA(ConcreteEJA): + +class CartesianProductEJA(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(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 = RealSymmetricEJA(3) - sage: J = DirectSumEJA(J1,J2) - sage: J.dimension() - 8 + 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 + 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 + + 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() + +----++----+----+----+ + | * || b0 | b1 | b2 | + +====++====+====+====+ + | b0 || b0 | 0 | 0 | + +----++----+----+----+ + | b1 || 0 | b1 | 0 | + +----++----+----+----+ + | b2 || 0 | 0 | b2 | + +----++----+----+----+ + sage: HadamardEJA(3).multiplication_table() + +----++----+----+----+ + | * || b0 | b1 | b2 | + +====++====+====+====+ + | b0 || b0 | 0 | 0 | + +----++----+----+----+ + | b1 || 0 | b1 | 0 | + +----++----+----+----+ + | b2 || 0 | 0 | b2 | + +----++----+----+----+ 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() - self._factors = (J1, J2) - basis = tuple( (a,b) for a in J1.basis() for b in J2.basis() ) + Element = FiniteDimensionalEJAElement - def jordan_product(x,y): - return (x[0]*y[0], x[1]*y[1]) - def inner_product(x,y): - return x[0].inner_product(y[0]) + x[1].inner_product(y[1]) + 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") - super().__init__(basis, jordan_product, inner_product) + associative = all( f.is_associative() for f in factors ) + + MS = self.matrix_space() + 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( + (factors[i](x[i])*factors[i](y[i])).to_matrix() + for i in range(m) + )) + + def inner_product(x, y): + return sum( + factors[i](x[i]).inner_product(factors[i](y[i])) + for i in range(m) + ) - self.rank.set_cache(J1.rank() + J2.rank()) + # 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().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 factors(self): + def matrix_space(self): r""" - Return the pair of this algebra's factors. + Return the space that our matrix basis lives in as a Cartesian + product. SETUP:: sage: from mjo.eja.eja_algebra import (HadamardEJA, - ....: JordanSpinEJA, - ....: DirectSumEJA) + ....: RealSymmetricEJA) EXAMPLES:: - sage: J1 = HadamardEJA(2, field=QQ) - sage: J2 = JordanSpinEJA(3, field=QQ) - sage: J = DirectSumEJA(J1,J2) - sage: J.factors() - (Euclidean Jordan algebra of dimension 2 over Rational Field, - Euclidean Jordan algebra of dimension 3 over Rational Field) - - """ - return self._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)) - - - -random_eja = ConcreteEJA.random_instance + 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 = 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 + + """ + 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 + def cartesian_embedding(self, i): + r""" + SETUP:: + + sage: from mjo.eja.eja_algebra import (random_eja, + ....: JordanSpinEJA, + ....: HadamardEJA, + ....: RealSymmetricEJA) + + 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 + + """ + 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 + +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