X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=ea62da8ca796031bd99382c5886bb1091d145b78;hb=bc02bf48592e22d034310cfffef8fb2a062c0a43;hp=12207b7c5a8e897738ab21a73361883cae03626f;hpb=5d646c586de50b571d2983b546a05899bf0c20c2;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 12207b7..b2891e5 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -1,103 +1,598 @@ """ -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` + * :class:`OctonionHermitianEJA` + +In addition to these, we provide two other example constructions, + + * :class:`JordanSpinEJA` + * :class:`HadamardEJA` + * :class:`AlbertEJA` + * :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. The Albert EJA is simply a special +case of the :class:`OctonionHermitianEJA` where the matrices are +three-by-three and the resulting space has dimension 27. And +last/least, the trivial EJA is exactly what you think it is; it could +also be obtained by constructing a dimension-zero instance of any of +the other algebras. 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. + +SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + +EXAMPLES:: + + sage: random_eja() + Euclidean Jordan algebra of dimension... """ -from itertools import izip, 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 from sage.misc.cachefunc import cached_method -from sage.misc.prandom import choice from sage.misc.table import table from sage.modules.free_module import FreeModule, VectorSpace -from sage.rings.all import (ZZ, QQ, RR, RLF, CLF, +from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF, PolynomialRing, QuadraticField) -from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement -from mjo.eja.eja_utils import _mat2vec - -class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): - # This is an ugly hack needed to prevent the category framework - # from implementing a coercion from our base ring (e.g. the - # rationals) into the algebra. First of all -- such a coercion is - # nonsense to begin with. But more importantly, it tries to do so - # in the category of rings, and since our algebras aren't - # associative they generally won't be rings. - _no_generic_basering_coercion = True +from mjo.eja.eja_element import FiniteDimensionalEJAElement +from mjo.eja.eja_operator import FiniteDimensionalEJAOperator +from mjo.eja.eja_utils import _all2list, _mat2vec + +class FiniteDimensionalEJA(CombinatorialFreeModule): + r""" + A finite-dimensional Euclidean Jordan algebra. + + INPUT: + + - ``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. + + - ``matrix_space`` -- the space that your matrix basis lives in, + or ``None`` (the default). So long as your basis does not have + length zero you can omit this. But in trivial algebras, it is + required. + + - ``field`` -- a subfield of the reals (default: ``AA``); the scalar + field for the algebra. + + - ``orthonormalize`` -- boolean (default: ``True``); whether or + not to orthonormalize the basis. Doing so is expensive and + generally rules out using the rationals as your ``field``, but + is required for spectral decompositions. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS: + + We should compute that an element subalgebra is associative even + if we circumvent the element method:: + + sage: set_random_seed() + sage: J = random_eja(field=QQ,orthonormalize=False) + sage: x = J.random_element() + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: basis = tuple(b.superalgebra_element() for b in A.basis()) + sage: J.subalgebra(basis, orthonormalize=False).is_associative() + True + """ + Element = FiniteDimensionalEJAElement def __init__(self, - field, - mult_table, - rank, - prefix='e', - category=None, - natural_basis=None, - check=True): + basis, + jordan_product, + inner_product, + field=AA, + matrix_space=None, + orthonormalize=True, + associative=None, + cartesian_product=False, + check_field=True, + check_axioms=True, + prefix="b"): + + n = len(basis) + + if check_field: + if not field.is_subring(RR): + # Note: this does return true for the real algebraic + # field, the rationals, and any quadratic field where + # we've specified a real embedding. + raise ValueError("scalar field is not real") + + if check_axioms: + # Check commutativity of the Jordan and inner-products. + # This has to be done before we build the multiplication + # and inner-product tables/matrices, because we take + # advantage of symmetry in the process. + if not all( jordan_product(bi,bj) == jordan_product(bj,bi) + for bi in basis + for bj in basis ): + raise ValueError("Jordan product is not commutative") + + if not all( inner_product(bi,bj) == inner_product(bj,bi) + for bi in basis + for bj in basis ): + raise ValueError("inner-product is not commutative") + + + category = MagmaticAlgebras(field).FiniteDimensional() + category = category.WithBasis().Unital().Commutative() + + if n <= 1: + # All zero- and one-dimensional algebras are just the real + # numbers with (some positive multiples of) the usual + # multiplication as its Jordan and inner-product. + associative = True + 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. + 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*b1 + 2*b2. + vector_basis = basis + + degree = 0 + if n > 0: + degree = len(_all2list(basis[0])) + + # Build an ambient space that fits our matrix basis when + # written out as "long vectors." + V = VectorSpace(field, degree) + + # The matrix that will hole the orthonormal -> unorthonormal + # coordinate transformation. + self._deortho_matrix = None + + if orthonormalize: + # 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(_all2list(b)) for b in basis ) + + from mjo.eja.eja_utils import gram_schmidt + basis = tuple(gram_schmidt(basis, inner_product)) + + # Save the (possibly orthonormalized) matrix basis for + # later, as well as the space that its elements live in. + # In most cases we can deduce the matrix space, but when + # n == 0 (that is, there are no basis elements) we cannot. + self._matrix_basis = basis + if matrix_space is None: + self._matrix_space = self._matrix_basis[0].parent() + else: + self._matrix_space = matrix_space + + # 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(_all2list(b)) for b in basis ) + W = V.span_of_basis( vector_basis, check=check_axioms) + + if orthonormalize: + # Now "W" is the vector space of our algebra coordinates. The + # variables "X1", "X2",... refer to the entries of vectors in + # W. Thus to convert back and forth between the orthonormal + # coordinates and the given ones, we need to stick the original + # basis in W. + U = V.span_of_basis( deortho_vector_basis, check=check_axioms) + self._deortho_matrix = matrix( U.coordinate_vector(q) + for q in vector_basis ) + + + # Now we actually compute the multiplication and inner-product + # tables/matrices using the possibly-orthonormalized basis. + self._inner_product_matrix = matrix.identity(field, n) + self._multiplication_table = [ [0 for j in range(i+1)] + for i in range(n) ] + + # Note: the Jordan and inner-products are defined in terms + # of the ambient basis. It's important that their arguments + # are in ambient coordinates as well. + for i in range(n): + for j in range(i+1): + # ortho basis w.r.t. ambient coords + q_i = basis[i] + q_j = basis[j] + + # 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(_all2list(elt))) + self._multiplication_table[i][j] = self.from_vector(elt) + + if not orthonormalize: + # If we're orthonormalizing the basis with respect + # to an inner-product, then the inner-product + # matrix with respect to the resulting basis is + # just going to be the identity. + ip = inner_product(q_i, q_j) + self._inner_product_matrix[i,j] = ip + self._inner_product_matrix[j,i] = ip + + self._inner_product_matrix._cache = {'hermitian': True} + self._inner_product_matrix.set_immutable() + + if check_axioms: + if not self._is_jordanian(): + raise ValueError("Jordan identity does not hold") + if not self._inner_product_is_associative(): + raise ValueError("inner product is not associative") + + + def _coerce_map_from_base_ring(self): + """ + Disable the map from the base ring into the algebra. + + Performing a nonsense conversion like this automatically + is counterpedagogical. The fallback is to try the usual + element constructor, which should also fail. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS:: + + sage: set_random_seed() + sage: J = random_eja() + sage: J(1) + Traceback (most recent call last): + ... + ValueError: not an element of this algebra + + """ + return None + + + 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: + return self._multiplication_table[i][j] + else: + return self._multiplication_table[j][i] + + def inner_product(self, x, y): """ + The inner product associated with this Euclidean Jordan algebra. + + Defaults to the trace inner product, but can be overridden by + subclasses if they are sure that the necessary properties are + satisfied. + SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja) + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, + ....: BilinearFormEJA) EXAMPLES: - By definition, Jordan multiplication commutes:: + Our inner product is "associative," which means the following for + a symmetric bilinear form:: sage: set_random_seed() sage: J = random_eja() - sage: x,y = J.random_elements(2) - sage: x*y == y*x + sage: x,y,z = J.random_elements(3) + sage: (x*y).inner_product(z) == y.inner_product(x*z) True TESTS: - The ``field`` we're given must be real:: + Ensure that this is the usual inner product for the algebras + over `R^n`:: - sage: JordanSpinEJA(2,QQbar) - Traceback (most recent call last): - ... - ValueError: field is not real + sage: set_random_seed() + sage: J = HadamardEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: actual = x.inner_product(y) + sage: expected = x.to_vector().inner_product(y.to_vector()) + sage: actual == expected + True + + Ensure that this is one-half of the trace inner-product in a + BilinearFormEJA that isn't just the reals (when ``n`` isn't + one). This is in Faraut and Koranyi, and also my "On the + symmetry..." paper:: + + sage: set_random_seed() + sage: J = BilinearFormEJA.random_instance() + sage: n = J.dimension() + sage: x = J.random_element() + sage: y = J.random_element() + sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) + True """ - if check: - if not field.is_subring(RR): - # Note: this does return true for the real algebraic - # field, and any quadratic field where we've specified - # a real embedding. - raise ValueError('field is not real') - - self._rank = rank - self._natural_basis = natural_basis - - if category is None: - category = MagmaticAlgebras(field).FiniteDimensional() - category = category.WithBasis().Unital() - - fda = super(FiniteDimensionalEuclideanJordanAlgebra, self) - fda.__init__(field, - range(len(mult_table)), - prefix=prefix, - category=category) - self.print_options(bracket='') - - # The multiplication table we're given is necessarily in terms - # of vectors, because we don't have an algebra yet for - # anything to be an element of. However, it's faster in the - # long run to have the multiplication table be in terms of - # algebra elements. We do this after calling the superclass - # constructor so that from_vector() knows what to do. - self._multiplication_table = [ map(lambda x: self.from_vector(x), ls) - for ls in mult_table ] + 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. + + This method should of course always return ``True``, unless + this algebra was constructed with ``check_axioms=False`` and + passed an invalid multiplication table. + """ + return all( x*y == y*x for x in self.gens() for y in self.gens() ) + + def _is_jordanian(self): + r""" + Whether or not this algebra's multiplication table respects the + Jordan identity `(x^{2})(xy) = x(x^{2}y)`. + + 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 + 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)) + == + (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j)) + for i in range(self.dimension()) + for j in range(self.dimension()) ) + + def _jordan_product_is_associative(self): + r""" + Return whether or not this algebra's Jordan product is + associative; that is, whether or not `x*(y*z) = (x*y)*z` + for all `x,y,x`. + This method should agree with :meth:`is_associative` unless + you lied about the value of the ``associative`` parameter + when you constructed the algebra. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (random_eja, + ....: RealSymmetricEJA, + ....: ComplexHermitianEJA, + ....: QuaternionHermitianEJA) + + EXAMPLES:: + + sage: J = RealSymmetricEJA(4, orthonormalize=False) + sage: J._jordan_product_is_associative() + False + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by() + sage: A._jordan_product_is_associative() + True + + :: + + sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) + sage: J._jordan_product_is_associative() + False + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: A._jordan_product_is_associative() + True + + :: + + sage: J = QuaternionHermitianEJA(2) + sage: J._jordan_product_is_associative() + False + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by() + sage: A._jordan_product_is_associative() + True + + TESTS: + + The values we've presupplied to the constructors agree with + the computation:: + + sage: set_random_seed() + sage: J = random_eja() + sage: J.is_associative() == J._jordan_product_is_associative() + True + + """ + R = self.base_ring() + + # Used to check whether or not something is zero. + epsilon = R.zero() + if not R.is_exact(): + # I don't know of any examples that make this magnitude + # necessary because I don't know how to make an + # associative algebra when the element subalgebra + # construction is unreliable (as it is over RDF; we can't + # find the degree of an element because we can't compute + # the rank of a matrix). But even multiplication of floats + # is non-associative, so *some* epsilon is needed... let's + # just take the one from _inner_product_is_associative? + epsilon = 1e-15 + + for i in range(self.dimension()): + for j in range(self.dimension()): + for k in range(self.dimension()): + x = self.monomial(i) + y = self.monomial(j) + z = self.monomial(k) + diff = (x*y)*z - x*(y*z) + + if diff.norm() > epsilon: + return False + + return True + + def _inner_product_is_associative(self): + r""" + Return whether or not this algebra's inner product `B` is + associative; that is, whether or not `B(xy,z) = B(x,yz)`. + + This method should of course always return ``True``, unless + 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. + 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) + diff = (x*y).inner_product(z) - x.inner_product(y*z) + + if diff.abs() > epsilon: + return False + + return True def _element_constructor_(self, elt): """ - Construct an element of this algebra from its natural + Construct an element of this algebra from its vector or matrix representation. This gets called only after the parent element _call_ method @@ -105,8 +600,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, - ....: RealCartesianProductEJA, + sage: from mjo.eja.eja_algebra import (random_eja, + ....: JordanSpinEJA, + ....: HadamardEJA, ....: RealSymmetricEJA) EXAMPLES: @@ -125,46 +621,86 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J(A) Traceback (most recent call last): ... - ArithmeticError: vector is not in free module + 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 natural representations are their usual - vector representations) 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 = RealCartesianProductEJA.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 """ - if elt == 0: - # The superclass implementation of random_element() - # needs to be able to coerce "0" into the algebra. - return self.zero() + msg = "not an element of this algebra" + 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 + # that the integer 3 belongs to the space of 2-by-2 matrices. + raise ValueError(msg) + + try: + # Try to convert a vector into a column-matrix... + elt = elt.column() + except (AttributeError, TypeError): + # and ignore failure, because we weren't really expecting + # a vector as an argument anyway. + pass - natural_basis = self.natural_basis() - basis_space = natural_basis[0].matrix_space() - if elt not in basis_space: - raise ValueError("not a naturally-represented algebra element") + if elt not in self.matrix_space(): + raise ValueError(msg) # Thanks for nothing! Matrix spaces aren't vector spaces in - # Sage, so we have to figure out its natural-basis coordinates + # Sage, so we have to figure out its matrix-basis coordinates # ourselves. We use the basis space's ring instead of the # element's ring because the basis space might be an algebraic # closure whereas the base ring of the 3-by-3 identity matrix # could be QQ instead of QQbar. - V = VectorSpace(basis_space.base_ring(), elt.nrows()*elt.ncols()) - W = V.span_of_basis( _mat2vec(s) for s in natural_basis ) - coords = W.coordinate_vector(_mat2vec(elt)) - return self.from_vector(coords) + # + # 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. + 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(V(elt)) + except ArithmeticError: # vector is not in free module + raise ValueError(msg) + return self.from_vector(coords) def _repr_(self): """ @@ -178,8 +714,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Ensure that it says what we think it says:: - sage: JordanSpinEJA(2, field=QQ) - Euclidean Jordan algebra of dimension 2 over Rational Field + sage: JordanSpinEJA(2, field=AA) + Euclidean Jordan algebra of dimension 2 over Algebraic Real Field sage: JordanSpinEJA(3, field=RDF) Euclidean Jordan algebra of dimension 3 over Real Double Field @@ -187,171 +723,14 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): fmt = "Euclidean Jordan algebra of dimension {} over {}" return fmt.format(self.dimension(), self.base_ring()) - def product_on_basis(self, i, j): - return self._multiplication_table[i][j] - - def _a_regular_element(self): - """ - Guess a regular element. Needed to compute the basis for our - characteristic polynomial coefficients. - - SETUP:: - - sage: from mjo.eja.eja_algebra import random_eja - - TESTS: - - Ensure that this hacky method succeeds for every algebra that we - know how to construct:: - - sage: set_random_seed() - sage: J = random_eja() - sage: J._a_regular_element().is_regular() - True - - """ - gs = self.gens() - z = self.sum( (i+1)*gs[i] for i in range(len(gs)) ) - if not z.is_regular(): - raise ValueError("don't know a regular element") - return z - - - @cached_method - def _charpoly_basis_space(self): - """ - Return the vector space spanned by the basis used in our - characteristic polynomial coefficients. This is used not only to - compute those coefficients, but also any time we need to - evaluate the coefficients (like when we compute the trace or - determinant). - """ - z = self._a_regular_element() - # Don't use the parent vector space directly here in case this - # happens to be a subalgebra. In that case, we would be e.g. - # two-dimensional but span_of_basis() would expect three - # coordinates. - V = VectorSpace(self.base_ring(), self.vector_space().dimension()) - basis = [ (z**k).to_vector() for k in range(self.rank()) ] - V1 = V.span_of_basis( basis ) - b = (V1.basis() + V1.complement().basis()) - return V.span_of_basis(b) - - - - @cached_method - def _charpoly_coeff(self, i): - """ - Return the coefficient polynomial "a_{i}" of this algebra's - general characteristic polynomial. - - Having this be a separate cached method lets us compute and - store the trace/determinant (a_{r-1} and a_{0} respectively) - separate from the entire characteristic polynomial. - """ - (A_of_x, x, xr, detA) = self._charpoly_matrix_system() - R = A_of_x.base_ring() - - if i == self.rank(): - return R.one() - if i > self.rank(): - # Guaranteed by theory - return R.zero() - - # Danger: the in-place modification is done for performance - # reasons (reconstructing a matrix with huge polynomial - # entries is slow), but I don't know how cached_method works, - # so it's highly possible that we're modifying some global - # list variable by reference, here. In other words, you - # probably shouldn't call this method twice on the same - # algebra, at the same time, in two threads - Ai_orig = A_of_x.column(i) - A_of_x.set_column(i,xr) - numerator = A_of_x.det() - A_of_x.set_column(i,Ai_orig) - - # We're relying on the theory here to ensure that each a_i is - # indeed back in R, and the added negative signs are to make - # the whole charpoly expression sum to zero. - return R(-numerator/detA) - - - @cached_method - def _charpoly_matrix_system(self): - """ - Compute the matrix whose entries A_ij are polynomials in - X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector - corresponding to `x^r` and the determinent of the matrix A = - [A_ij]. In other words, all of the fixed (cachable) data needed - to compute the coefficients of the characteristic polynomial. - """ - r = self.rank() - n = self.dimension() - - # Turn my vector space into a module so that "vectors" can - # have multivatiate polynomial entries. - names = tuple('X' + str(i) for i in range(1,n+1)) - R = PolynomialRing(self.base_ring(), names) - - # Using change_ring() on the parent's vector space doesn't work - # here because, in a subalgebra, that vector space has a basis - # and change_ring() tries to bring the basis along with it. And - # that doesn't work unless the new ring is a PID, which it usually - # won't be. - V = FreeModule(R,n) - - # Now let x = (X1,X2,...,Xn) be the vector whose entries are - # indeterminates... - x = V(names) - - # And figure out the "left multiplication by x" matrix in - # that setting. - lmbx_cols = [] - monomial_matrices = [ self.monomial(i).operator().matrix() - for i in range(n) ] # don't recompute these! - for k in range(n): - ek = self.monomial(k).to_vector() - lmbx_cols.append( - sum( x[i]*(monomial_matrices[i]*ek) - for i in range(n) ) ) - Lx = matrix.column(R, lmbx_cols) - - # Now we can compute powers of x "symbolically" - x_powers = [self.one().to_vector(), x] - for d in range(2, r+1): - x_powers.append( Lx*(x_powers[-1]) ) - - idmat = matrix.identity(R, n) - - W = self._charpoly_basis_space() - W = W.change_ring(R.fraction_field()) - - # Starting with the standard coordinates x = (X1,X2,...,Xn) - # and then converting the entries to W-coordinates allows us - # to pass in the standard coordinates to the charpoly and get - # back the right answer. Specifically, with x = (X1,X2,...,Xn), - # we have - # - # W.coordinates(x^2) eval'd at (standard z-coords) - # = - # W-coords of (z^2) - # = - # W-coords of (standard coords of x^2 eval'd at std-coords of z) - # - # We want the middle equivalent thing in our matrix, but use - # the first equivalent thing instead so that we can pass in - # standard coordinates. - x_powers = [ W.coordinate_vector(xp) for xp in x_powers ] - l2 = [idmat.column(k-1) for k in range(r+1, n+1)] - A_of_x = matrix.column(R, n, (x_powers[:r] + l2)) - return (A_of_x, x, x_powers[r], A_of_x.det()) - @cached_method - def characteristic_polynomial(self): + def characteristic_polynomial_of(self): """ - Return a characteristic polynomial that works for all elements - of this algebra. + Return the algebra's "characteristic polynomial of" function, + which is itself a multivariate polynomial that, when evaluated + at the coordinates of some algebra element, returns that + element's characteristic polynomial. The resulting polynomial has `n+1` variables, where `n` is the dimension of this algebra. The first `n` variables correspond to @@ -371,7 +750,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Alizadeh, Example 11.11:: sage: J = JordanSpinEJA(3) - sage: p = J.characteristic_polynomial(); p + sage: p = J.characteristic_polynomial_of(); p X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2 sage: xvec = J.one().to_vector() sage: p(*xvec) @@ -384,28 +763,51 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): any argument:: sage: J = TrivialEJA() - sage: J.characteristic_polynomial() + sage: J.characteristic_polynomial_of() 1 """ r = self.rank() n = self.dimension() - # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n. - a = [ self._charpoly_coeff(i) for i in range(r+1) ] + # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1). + a = self._charpoly_coefficients() # We go to a bit of trouble here to reorder the # indeterminates, so that it's easier to evaluate the # characteristic polynomial at x's coordinates and get back # something in terms of t, which is what we want. - R = a[0].parent() S = PolynomialRing(self.base_ring(),'t') t = S.gen(0) - S = PolynomialRing(S, R.variable_names()) - t = S(t) + if r > 0: + R = a[0].parent() + S = PolynomialRing(S, R.variable_names()) + t = S(t) + + return (t**r + sum( a[k]*(t**k) for k in range(r) )) + + def coordinate_polynomial_ring(self): + r""" + The multivariate polynomial ring in which this algebra's + :meth:`characteristic_polynomial_of` lives. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: RealSymmetricEJA) + + EXAMPLES:: - return sum( a[k]*(t**k) for k in xrange(len(a)) ) + sage: J = HadamardEJA(2) + sage: J.coordinate_polynomial_ring() + Multivariate Polynomial Ring in X1, X2... + sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False) + sage: J.coordinate_polynomial_ring() + Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6... + """ + var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) ) + return PolynomialRing(self.base_ring(), var_names) def inner_product(self, x, y): """ @@ -417,7 +819,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): SETUP:: - sage: from mjo.eja.eja_algebra import random_eja + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, + ....: BilinearFormEJA) EXAMPLES: @@ -430,10 +834,34 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: (x*y).inner_product(z) == y.inner_product(x*z) True + TESTS: + + Ensure that this is the usual inner product for the algebras + over `R^n`:: + + sage: set_random_seed() + sage: J = HadamardEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: actual = x.inner_product(y) + sage: expected = x.to_vector().inner_product(y.to_vector()) + sage: actual == expected + True + + Ensure that this is one-half of the trace inner-product in a + BilinearFormEJA that isn't just the reals (when ``n`` isn't + one). This is in Faraut and Koranyi, and also my "On the + symmetry..." paper:: + + sage: set_random_seed() + sage: J = BilinearFormEJA.random_instance() + sage: n = J.dimension() + sage: x = J.random_element() + sage: y = J.random_element() + sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) + True """ - X = x.natural_representation() - Y = y.natural_representation() - return self.natural_inner_product(X,Y) + B = self._inner_product_matrix + return (B*x.to_vector()).inner_product(y.to_vector()) def is_trivial(self): @@ -477,38 +905,58 @@ class FiniteDimensionalEuclideanJordanAlgebra(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 | +----++----+----+----+----+ """ - M = list(self._multiplication_table) # copy - for i in xrange(len(M)): - # M had better be "square" - M[i] = [self.monomial(i)] + M[i] - M = [["*"] + list(self.gens())] + M + n = self.dimension() + # Prepend the header row. + M = [["*"] + list(self.gens())] + + # And to each subsequent row, prepend an entry that belongs to + # the left-side "header column." + 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) - def natural_basis(self): + def matrix_basis(self): """ - Return a more-natural representation of this algebra's basis. + Return an (often more natural) representation of this algebras + basis as an ordered tuple of matrices. - Every finite-dimensional Euclidean Jordan Algebra is a direct - sum of five simple algebras, four of which comprise Hermitian - matrices. This method returns the original "natural" basis - for our underlying vector space. (Typically, the natural basis - is used to construct the multiplication table in the first place.) + Every finite-dimensional Euclidean Jordan Algebra is a, up to + Jordan isomorphism, a direct sum of five simple + algebras---four of which comprise Hermitian matrices. And the + last type of algebra can of course be thought of as `n`-by-`1` + column matrices (ambiguusly called column vectors) to avoid + special cases. As a result, matrices (and column vectors) are + a natural representation format for Euclidean Jordan algebra + elements. - Note that this will always return a matrix. The standard basis - in `R^n` will be returned as `n`-by-`1` column matrices. + But, when we construct an algebra from a basis of matrices, + those matrix representations are lost in favor of coordinate + vectors *with respect to* that basis. We could eventually + convert back if we tried hard enough, but having the original + representations handy is valuable enough that we simply store + them and return them from this method. + + 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:`CartesianProductEJA` can be displayed + nicely, without having to have special classes for direct sums + one of whose components was a matrix algebra. SETUP:: @@ -519,55 +967,81 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = RealSymmetricEJA(2) sage: J.basis() - Finite family {0: e0, 1: e1, 2: e2} - sage: J.natural_basis() + Finite family {0: b0, 1: b1, 2: b2} + sage: J.matrix_basis() ( - [1 0] [ 0 1/2*sqrt2] [0 0] - [0 0], [1/2*sqrt2 0], [0 1] + [1 0] [ 0 0.7071067811865475?] [0 0] + [0 0], [0.7071067811865475? 0], [0 1] ) :: sage: J = JordanSpinEJA(2) sage: J.basis() - Finite family {0: e0, 1: e1} - sage: J.natural_basis() + Finite family {0: b0, 1: b1} + sage: J.matrix_basis() ( [1] [0] [0], [1] ) - """ - if self._natural_basis is None: - M = self.natural_basis_space() - return tuple( M(b.to_vector()) for b in self.basis() ) - else: - return self._natural_basis + return self._matrix_basis - def natural_basis_space(self): - """ - Return the matrix space in which this algebra's natural basis - elements live. + def matrix_space(self): """ - if self._natural_basis is None or len(self._natural_basis) == 0: - return MatrixSpace(self.base_ring(), self.dimension(), 1) - else: - return self._natural_basis[0].matrix_space() + Return the matrix space in which this algebra's elements live, if + we think of them as matrices (including column vectors of the + appropriate size). + "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. For algebras of + matrices, this returns the space in which their + real embeddings live. - @staticmethod - def natural_inner_product(X,Y): - """ - Compute the inner product of two naturally-represented elements. + 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() + Module of 2 by 2 matrices with entries in Algebraic Field over + the scalar ring Rational Field + sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False) + sage: J.matrix_space() + Module of 1 by 1 matrices with entries in Quaternion + Algebra (-1, -1) with base ring Rational Field over + the scalar ring Rational Field - For example in the real symmetric matrix EJA, this will compute - the trace inner-product of two n-by-n symmetric matrices. The - default should work for the real cartesian product EJA, the - Jordan spin EJA, and the real symmetric matrices. The others - will have to be overridden. """ - return (X.conjugate_transpose()*Y).trace() + return self._matrix_space @cached_method @@ -577,26 +1051,60 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): SETUP:: - sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA, + sage: from mjo.eja.eja_algebra import (HadamardEJA, ....: random_eja) - EXAMPLES:: + EXAMPLES: - sage: J = RealCartesianProductEJA(5) + We can compute unit element in the Hadamard EJA:: + + sage: J = HadamardEJA(5) + sage: J.one() + 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() + c0 + sage: A.one().superalgebra_element() + b0 + b1 + b2 + b3 + b4 TESTS: - The identity element acts like the identity:: + The identity element acts like the identity, regardless of + whether or not we orthonormalize:: sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: J.one()*x == x and x*J.one() == x True + sage: A = x.subalgebra_generated_by() + sage: y = A.random_element() + sage: A.one()*y == y and y*A.one() == y + True + + :: + + sage: set_random_seed() + sage: J = random_eja(field=QQ, orthonormalize=False) + sage: x = J.random_element() + sage: J.one()*x == x and x*J.one() == x + True + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: y = A.random_element() + sage: A.one()*y == y and y*A.one() == y + True - The matrix of the unit element's operator is the identity:: + The matrix of the unit element's operator is the identity, + regardless of the base field and whether or not we + orthonormalize:: sage: set_random_seed() sage: J = random_eja() @@ -604,37 +1112,272 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: expected = matrix.identity(J.base_ring(), J.dimension()) sage: actual == expected True + sage: x = J.random_element() + sage: A = x.subalgebra_generated_by() + sage: actual = A.one().operator().matrix() + sage: expected = matrix.identity(A.base_ring(), A.dimension()) + sage: actual == expected + True - """ - # We can brute-force compute the matrices of the operators - # that correspond to the basis elements of this algebra. - # If some linear combination of those basis elements is the - # algebra identity, then the same linear combination of - # their matrices has to be the identity matrix. - # - # Of course, matrices aren't vectors in sage, so we have to + :: + + sage: set_random_seed() + sage: J = random_eja(field=QQ, orthonormalize=False) + sage: actual = J.one().operator().matrix() + sage: expected = matrix.identity(J.base_ring(), J.dimension()) + sage: actual == expected + True + sage: x = J.random_element() + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: actual = A.one().operator().matrix() + sage: expected = matrix.identity(A.base_ring(), A.dimension()) + sage: actual == expected + True + + Ensure that the cached unit element (often precomputed by + hand) agrees with the computed one:: + + sage: set_random_seed() + sage: J = random_eja() + sage: cached = J.one() + sage: J.one.clear_cache() + sage: J.one() == cached + True + + :: + + sage: set_random_seed() + sage: J = random_eja(field=QQ, orthonormalize=False) + sage: cached = J.one() + sage: J.one.clear_cache() + sage: J.one() == cached + True + + """ + # We can brute-force compute the matrices of the operators + # that correspond to the basis elements of this algebra. + # If some linear combination of those basis elements is the + # algebra identity, then the same linear combination of + # their matrices has to be the identity matrix. + # + # Of course, matrices aren't vectors in sage, so we have to # appeal to the "long vectors" isometry. oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ] - # Now we use basis linear algebra to find the coefficients, + # Now we use basic linear algebra to find the coefficients, # of the matrices-as-vectors-linear-combination, which should # work for the original algebra basis too. - A = matrix.column(self.base_ring(), oper_vecs) + A = matrix(self.base_ring(), oper_vecs) # We used the isometry on the left-hand side already, but we # still need to do it for the right-hand side. Recall that we # wanted something that summed to the identity matrix. b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) ) - # Now if there's an identity element in the algebra, this should work. - coeffs = A.solve_right(b) - return self.linear_combination(zip(self.gens(), coeffs)) + # Now if there's an identity element in the algebra, this + # should work. We solve on the left to avoid having to + # transpose the matrix "A". + return self.from_vector(A.solve_left(b)) + + + def peirce_decomposition(self, c): + """ + The Peirce decomposition of this algebra relative to the + idempotent ``c``. + + In the future, this can be extended to a complete system of + orthogonal idempotents. + + INPUT: + + - ``c`` -- an idempotent of this algebra. + + OUTPUT: + + A triple (J0, J5, J1) containing two subalgebras and one subspace + of this algebra, + + - ``J0`` -- the algebra on the eigenspace of ``c.operator()`` + corresponding to the eigenvalue zero. + + - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()`` + corresponding to the eigenvalue one-half. + + - ``J1`` -- the algebra on the eigenspace of ``c.operator()`` + corresponding to the eigenvalue one. + + These are the only possible eigenspaces for that operator, and this + algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are + orthogonal, and are subalgebras of this algebra with the appropriate + restrictions. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA + + EXAMPLES: + + The canonical example comes from the symmetric matrices, which + decompose into diagonal and off-diagonal parts:: + + sage: J = RealSymmetricEJA(3) + sage: C = matrix(QQ, [ [1,0,0], + ....: [0,1,0], + ....: [0,0,0] ]) + sage: c = J(C) + sage: J0,J5,J1 = J.peirce_decomposition(c) + sage: J0 + Euclidean Jordan algebra of dimension 1... + sage: J5 + Vector space of degree 6 and dimension 2... + sage: J1 + Euclidean Jordan algebra of dimension 3... + sage: J0.one().to_matrix() + [0 0 0] + [0 0 0] + [0 0 1] + sage: orig_df = AA.options.display_format + sage: AA.options.display_format = 'radical' + sage: J.from_vector(J5.basis()[0]).to_matrix() + [ 0 0 1/2*sqrt(2)] + [ 0 0 0] + [1/2*sqrt(2) 0 0] + sage: J.from_vector(J5.basis()[1]).to_matrix() + [ 0 0 0] + [ 0 0 1/2*sqrt(2)] + [ 0 1/2*sqrt(2) 0] + sage: AA.options.display_format = orig_df + sage: J1.one().to_matrix() + [1 0 0] + [0 1 0] + [0 0 0] + + TESTS: + + Every algebra decomposes trivially with respect to its identity + element:: + + sage: set_random_seed() + sage: J = random_eja() + sage: J0,J5,J1 = J.peirce_decomposition(J.one()) + sage: J0.dimension() == 0 and J5.dimension() == 0 + True + sage: J1.superalgebra() == J and J1.dimension() == J.dimension() + True + + The decomposition is into eigenspaces, and its components are + therefore necessarily orthogonal. Moreover, the identity + elements in the two subalgebras are the projections onto their + respective subspaces of the superalgebra's identity element:: + + sage: set_random_seed() + sage: J = random_eja() + sage: x = J.random_element() + sage: if not J.is_trivial(): + ....: while x.is_nilpotent(): + ....: x = J.random_element() + sage: c = x.subalgebra_idempotent() + sage: J0,J5,J1 = J.peirce_decomposition(c) + sage: ipsum = 0 + sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()): + ....: w = w.superalgebra_element() + ....: y = J.from_vector(y) + ....: z = z.superalgebra_element() + ....: ipsum += w.inner_product(y).abs() + ....: ipsum += w.inner_product(z).abs() + ....: ipsum += y.inner_product(z).abs() + sage: ipsum + 0 + sage: J1(c) == J1.one() + True + sage: J0(J.one() - c) == J0.one() + True + + """ + if not c.is_idempotent(): + raise ValueError("element is not idempotent: %s" % c) + + # 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 = self.subalgebra(()) + J0 = trivial # eigenvalue zero + J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half + J1 = trivial # eigenvalue one + + for (eigval, eigspace) in c.operator().matrix().right_eigenspaces(): + if eigval == ~(self.base_ring()(2)): + J5 = eigspace + else: + gens = tuple( self.from_vector(b) for b in eigspace.basis() ) + subalg = self.subalgebra(gens, check_axioms=False) + if eigval == 0: + J0 = subalg + elif eigval == 1: + J1 = subalg + else: + raise ValueError("unexpected eigenvalue: %s" % eigval) + + return (J0, J5, J1) + + + def random_element(self, thorough=False): + r""" + Return a random element of this algebra. + + Our algebra superclass method only returns a linear + combination of at most two basis elements. We instead + want the vector space "random element" method that + returns a more diverse selection. + + INPUT: + + - ``thorough`` -- (boolean; default False) whether or not we + should generate irrational coefficients for the random + element when our base ring is irrational; this slows the + algebra operations to a crawl, but any truly random method + should include them + + """ + # For a general base ring... maybe we can trust this to do the + # right thing? Unlikely, but. + V = self.vector_space() + v = V.random_element() + + if self.base_ring() is AA: + # The "random element" method of the algebraic reals is + # stupid at the moment, and only returns integers between + # -2 and 2, inclusive: + # + # https://trac.sagemath.org/ticket/30875 + # + # Instead, we implement our own "random vector" method, + # and then coerce that into the algebra. We use the vector + # space degree here instead of the dimension because a + # subalgebra could (for example) be spanned by only two + # vectors, each with five coordinates. We need to + # generate all five coordinates. + if thorough: + v *= QQbar.random_element().real() + else: + v *= QQ.random_element() + return self.from_vector(V.coordinate_vector(v)) - def random_elements(self, count): + def random_elements(self, count, thorough=False): """ Return ``count`` random elements as a tuple. + INPUT: + + - ``thorough`` -- (boolean; default False) whether or not we + should generate irrational coefficients for the random + elements when our base ring is irrational; this slows the + algebra operations to a crawl, but any truly random method + should include them + SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA @@ -649,23 +1392,91 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): True """ - return tuple( self.random_element() for idx in xrange(count) ) + return tuple( self.random_element(thorough) + for idx in range(count) ) - def rank(self): + @cached_method + def _charpoly_coefficients(self): + 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 + """ - Return the rank of this EJA. + n = self.dimension() + R = self.coordinate_polynomial_ring() + vars = R.gens() + F = R.fraction_field() + + 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] + for k in range(n) ) + + L_x = matrix(F, n, n, L_x_i_j) + + r = None + if self.rank.is_in_cache(): + r = self.rank() + # There's no need to pad the system with redundant + # columns if we *know* they'll be redundant. + n = r + + # Compute an extra power in case the rank is equal to + # the dimension (otherwise, we would stop at x^(r-1)). + x_powers = [ (L_x**k)*self.one().to_vector() + for k in range(n+1) ] + A = matrix.column(F, x_powers[:n]) + AE = A.extended_echelon_form() + E = AE[:,n:] + A_rref = AE[:,:n] + if r is None: + r = A_rref.rank() + b = x_powers[r] + + # 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. 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) - ALGORITHM: + @cached_method + def rank(self): + r""" + Return the rank of this EJA. - The author knows of no algorithm to compute the rank of an EJA - where only the multiplication table is known. In lieu of one, we - require the rank to be specified when the algebra is created, - and simply pass along that number here. + This is a cached method because we know the rank a priori for + all of the algebras we can construct. Thus we can avoid the + expensive ``_charpoly_coefficients()`` call unless we truly + need to compute the whole characteristic polynomial. SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: JordanSpinEJA, ....: RealSymmetricEJA, ....: ComplexHermitianEJA, ....: QuaternionHermitianEJA, @@ -706,8 +1517,26 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: r > 0 or (r == 0 and J.is_trivial()) True + Ensure that computing the rank actually works, since the ranks + of all simple algebras are known and will be cached by default:: + + sage: set_random_seed() # long time + sage: J = random_eja() # long time + sage: cached = J.rank() # long time + sage: J.rank.clear_cache() # long time + sage: J.rank() == cached # long time + True + + """ + return len(self._charpoly_coefficients()) + + + def subalgebra(self, basis, **kwargs): + r""" + Create a subalgebra of this algebra from the given basis. """ - return self._rank + from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra + return FiniteDimensionalEJASubalgebra(self, basis, **kwargs) def vector_space(self): @@ -728,335 +1557,353 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return self.zero().to_vector().parent().ambient_vector_space() - Element = FiniteDimensionalEuclideanJordanAlgebraElement +class RationalBasisEJA(FiniteDimensionalEJA): + r""" + Algebras whose supplied basis elements have all rational entries. -class KnownRankEJA(object): - """ - A class for algebras that we actually know we can construct. The - main issue is that, for most of our methods to make sense, we need - to know the rank of our algebra. Thus we can't simply generate a - "random" algebra, or even check that a given basis and product - satisfy the axioms; because even if everything looks OK, we wouldn't - know the rank we need to actuallty build the thing. + SETUP:: - Not really a subclass of FDEJA because doing that causes method - resolution errors, e.g. + sage: from mjo.eja.eja_algebra import BilinearFormEJA - TypeError: Error when calling the metaclass bases - Cannot create a consistent method resolution - order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra, - KnownRankEJA + EXAMPLES: - """ - @staticmethod - def _max_test_case_size(): - """ - Return an integer "size" that is an upper bound on the size of - this algebra when it is used in a random test - case. Unfortunately, the term "size" is quite vague -- when - dealing with `R^n` under either the Hadamard or Jordan spin - product, the "size" refers to the dimension `n`. When dealing - with a matrix algebra (real symmetric or complex/quaternion - Hermitian), it refers to the size of the matrix, which is - far less than the dimension of the underlying vector space. + The supplied basis is orthonormalized by default:: - We default to five in this class, which is safe in `R^n`. The - matrix algebra subclasses (or any class where the "size" is - interpreted to be far less than the dimension) should override - with a smaller number. - """ - return 5 + sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]]) + sage: J = BilinearFormEJA(B) + sage: J.matrix_basis() + ( + [1] [ 0] [ 0] + [0] [1/5] [32/5] + [0], [ 0], [ 5] + ) - @classmethod - def random_instance(cls, field=QQ, **kwargs): - """ - Return a random instance of this type of algebra. + """ + def __init__(self, + basis, + jordan_product, + inner_product, + field=AA, + check_field=True, + **kwargs): + + if check_field: + # Abuse the check_field parameter to check that the entries of + # out basis (in ambient coordinates) are in the field QQ. + # Use _all2list to get the vector coordinates of octonion + # entries and not the octonions themselves (which are not + # rational). + if not all( all(b_i in QQ for b_i in _all2list(b)) + 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 + # one is already rational. + # + # Note: the same Jordan and inner-products work here, + # because they are necessarily defined with respect to + # ambient coordinates and not any particular basis. + self._rational_algebra = FiniteDimensionalEJA( + basis, + jordan_product, + inner_product, + field=QQ, + matrix_space=self.matrix_space(), + associative=self.is_associative(), + orthonormalize=False, + check_field=False, + check_axioms=False) - Beware, this will crash for "most instances" because the - constructor below looks wrong. - """ - if cls is TrivialEJA: - # The TrivialEJA class doesn't take an "n" argument because - # there's only one. - return cls(field) + @cached_method + def _charpoly_coefficients(self): + r""" + SETUP:: - n = ZZ.random_element(cls._max_test_case_size()) + 1 - return cls(n, field, **kwargs) + sage: from mjo.eja.eja_algebra import (BilinearFormEJA, + ....: JordanSpinEJA) + EXAMPLES: -class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra, - KnownRankEJA): - """ - Return the Euclidean Jordan Algebra corresponding to the set - `R^n` under the Hadamard product. + The base ring of the resulting polynomial coefficients is what + it should be, and not the rationals (unless the algebra was + already over the rationals):: - Note: this is nothing more than the Cartesian product of ``n`` - copies of the spin algebra. Once Cartesian product algebras - are implemented, this can go. + sage: J = JordanSpinEJA(3) + sage: J._charpoly_coefficients() + (X1^2 - X2^2 - X3^2, -2*X1) + sage: a0 = J._charpoly_coefficients()[0] + sage: J.base_ring() + Algebraic Real Field + sage: a0.base_ring() + Algebraic Real Field + + """ + if self._rational_algebra is None: + # There's no need to construct *another* algebra over the + # rationals if this one is already over the + # rationals. Likewise, if we never orthonormalized our + # basis, we might as well just use the given one. + return super()._charpoly_coefficients() + + # Do the computation over the rationals. The answer will be + # the same, because all we've done is a change of basis. + # Then, change back from QQ to our real base ring + a = ( a_i.change_ring(self.base_ring()) + for a_i in self._rational_algebra._charpoly_coefficients() ) + + if self._deortho_matrix is None: + # This can happen if our base ring was, say, AA and we + # chose not to (or didn't need to) orthonormalize. It's + # still faster to do the computations over QQ even if + # the numbers in the boxes stay the same. + return tuple(a) + + # Otherwise, convert the coordinate variables back to the + # deorthonormalized ones. + R = self.coordinate_polynomial_ring() + from sage.modules.free_module_element import vector + X = vector(R, R.gens()) + BX = self._deortho_matrix*X + + subs_dict = { X[i]: BX[i] for i in range(len(X)) } + return tuple( a_i.subs(subs_dict) for a_i in a ) + +class ConcreteEJA(FiniteDimensionalEJA): + r""" + A class for the Euclidean Jordan algebras that we know by name. + + These are the Jordan algebras whose basis, multiplication table, + rank, and so on are known a priori. More to the point, they are + the Euclidean Jordan algebras for which we are able to conjure up + a "random instance." SETUP:: - sage: from mjo.eja.eja_algebra import RealCartesianProductEJA - - EXAMPLES: - - This multiplication table can be verified by hand:: - - sage: J = RealCartesianProductEJA(3) - sage: e0,e1,e2 = J.gens() - sage: e0*e0 - e0 - sage: e0*e1 - 0 - sage: e0*e2 - 0 - sage: e1*e1 - e1 - sage: e1*e2 - 0 - sage: e2*e2 - e2 + sage: from mjo.eja.eja_algebra import ConcreteEJA TESTS: - We can change the generator prefix:: + Our basis is normalized with respect to the algebra's inner + product, unless we specify otherwise:: - sage: RealCartesianProductEJA(3, prefix='r').gens() - (r0, r1, r2) + sage: set_random_seed() + sage: J = ConcreteEJA.random_instance() + sage: all( b.norm() == 1 for b in J.gens() ) + True - """ - def __init__(self, n, field=QQ, **kwargs): - V = VectorSpace(field, n) - mult_table = [ [ V.gen(i)*(i == j) for j in xrange(n) ] - for i in xrange(n) ] + Since our basis is orthonormal with respect to the algebra's inner + product, and since we know that this algebra is an EJA, any + left-multiplication operator's matrix will be symmetric because + natural->EJA basis representation is an isometry and within the + EJA the operator is self-adjoint by the Jordan axiom:: - fdeja = super(RealCartesianProductEJA, self) - return fdeja.__init__(field, mult_table, rank=n, **kwargs) + sage: set_random_seed() + sage: J = ConcreteEJA.random_instance() + sage: x = J.random_element() + sage: x.operator().is_self_adjoint() + True + """ - def inner_product(self, x, y): + @staticmethod + def _max_random_instance_size(): """ - Faster to reimplement than to use natural representations. + Return an integer "size" that is an upper bound on the size of + this algebra when it is used in a random test + case. Unfortunately, the term "size" is ambiguous -- when + dealing with `R^n` under either the Hadamard or Jordan spin + product, the "size" refers to the dimension `n`. When dealing + with a matrix algebra (real symmetric or complex/quaternion + Hermitian), it refers to the size of the matrix, which is far + less than the dimension of the underlying vector space. - SETUP:: + This method must be implemented in each subclass. + """ + raise NotImplementedError - sage: from mjo.eja.eja_algebra import RealCartesianProductEJA + @classmethod + def random_instance(cls, *args, **kwargs): + """ + Return a random instance of this type of algebra. - TESTS: + This method should be implemented in each subclass. + """ + from sage.misc.prandom import choice + eja_class = choice(cls.__subclasses__()) - Ensure that this is the usual inner product for the algebras - over `R^n`:: + # These all bubble up to the RationalBasisEJA superclass + # constructor, so any (kw)args valid there are also valid + # here. + return eja_class.random_instance(*args, **kwargs) - sage: set_random_seed() - sage: J = RealCartesianProductEJA.random_instance() - sage: x,y = J.random_elements(2) - sage: X = x.natural_representation() - sage: Y = y.natural_representation() - sage: x.inner_product(y) == J.natural_inner_product(X,Y) - True +class MatrixEJA(FiniteDimensionalEJA): + @staticmethod + def _denormalized_basis(A): """ - return x.to_vector().inner_product(y.to_vector()) - + Returns a basis for the space of complex Hermitian n-by-n matrices. -def random_eja(field=QQ, nontrivial=False): - """ - Return a "random" finite-dimensional Euclidean Jordan Algebra. + Why do we embed these? Basically, because all of numerical linear + algebra assumes that you're working with vectors consisting of `n` + entries from a field and scalars from the same field. There's no way + to tell SageMath that (for example) the vectors contain complex + numbers, while the scalar field is real. - SETUP:: + SETUP:: - sage: from mjo.eja.eja_algebra import random_eja + sage: from mjo.hurwitz import (ComplexMatrixAlgebra, + ....: QuaternionMatrixAlgebra, + ....: OctonionMatrixAlgebra) + sage: from mjo.eja.eja_algebra import MatrixEJA - TESTS:: + TESTS:: - sage: random_eja() - Euclidean Jordan algebra of dimension... + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: A = MatrixSpace(QQ, n) + sage: B = MatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B) + True - """ - eja_classes = KnownRankEJA.__subclasses__() - if nontrivial: - eja_classes.remove(TrivialEJA) - classname = choice(eja_classes) - return classname.random_instance(field=field) + :: + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: A = ComplexMatrixAlgebra(n, scalars=QQ) + sage: B = MatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B) + True + :: + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: A = QuaternionMatrixAlgebra(n, scalars=QQ) + sage: B = MatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B ) + True + :: + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: A = OctonionMatrixAlgebra(n, scalars=QQ) + sage: B = MatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B ) + True -class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): - @staticmethod - def _max_test_case_size(): - # Play it safe, since this will be squared and the underlying - # field can have dimension 4 (quaternions) too. - return 2 - - def __init__(self, field, basis, rank, normalize_basis=True, **kwargs): - """ - Compared to the superclass constructor, we take a basis instead of - a multiplication table because the latter can be computed in terms - of the former when the product is known (like it is here). - """ - # Used in this class's fast _charpoly_coeff() override. - self._basis_normalizers = None - - # We're going to loop through this a few times, so now's a good - # time to ensure that it isn't a generator expression. - basis = tuple(basis) - - if rank > 1 and normalize_basis: - # We'll need sqrt(2) to normalize the basis, and this - # winds up in the multiplication table, so the whole - # algebra needs to be over the field extension. - R = PolynomialRing(field, 'z') - z = R.gen() - p = z**2 - 2 - if p.is_irreducible(): - field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt()) - basis = tuple( s.change_ring(field) for s in basis ) - self._basis_normalizers = tuple( - ~(self.natural_inner_product(s,s).sqrt()) for s in basis ) - basis = tuple(s*c for (s,c) in izip(basis,self._basis_normalizers)) - - Qs = self.multiplication_table_from_matrix_basis(basis) - - fdeja = super(MatrixEuclideanJordanAlgebra, self) - return fdeja.__init__(field, - Qs, - rank=rank, - natural_basis=basis, - **kwargs) + """ + # These work for real MatrixSpace, whose monomials only have + # two coordinates (because the last one would always be "1"). + es = A.base_ring().gens() + gen = lambda A,m: A.monomial(m[:2]) + if hasattr(A, 'entry_algebra_gens'): + # We've got a MatrixAlgebra, and its monomials will have + # three coordinates. + es = A.entry_algebra_gens() + gen = lambda A,m: A.monomial(m) - @cached_method - def _charpoly_coeff(self, i): - """ - Override the parent method with something that tries to compute - over a faster (non-extension) field. - """ - if self._basis_normalizers is None: - # We didn't normalize, so assume that the basis we started - # with had entries in a nice field. - return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i) - else: - basis = ( (b/n) for (b,n) in izip(self.natural_basis(), - self._basis_normalizers) ) - - # Do this over the rationals and convert back at the end. - J = MatrixEuclideanJordanAlgebra(QQ, - basis, - self.rank(), - normalize_basis=False) - (_,x,_,_) = J._charpoly_matrix_system() - p = J._charpoly_coeff(i) - # p might be missing some vars, have to substitute "optionally" - pairs = izip(x.base_ring().gens(), self._basis_normalizers) - substitutions = { v: v*c for (v,c) in pairs } - result = p.subs(substitutions) - - # The result of "subs" can be either a coefficient-ring - # element or a polynomial. Gotta handle both cases. - if result in QQ: - return self.base_ring()(result) - else: - return result.change_ring(self.base_ring()) + basis = [] + for i in range(A.nrows()): + for j in range(i+1): + if i == j: + E_ii = gen(A, (i,j,es[0])) + basis.append(E_ii) + else: + for e in es: + E_ij = gen(A, (i,j,e)) + E_ij += E_ij.conjugate_transpose() + basis.append(E_ij) + return tuple( basis ) @staticmethod - def multiplication_table_from_matrix_basis(basis): - """ - At least three of the five simple Euclidean Jordan algebras have the - symmetric multiplication (A,B) |-> (AB + BA)/2, where the - multiplication on the right is matrix multiplication. Given a basis - for the underlying matrix space, this function returns a - multiplication table (obtained by looping through the basis - elements) for an algebra of those matrices. - """ - # In S^2, for example, we nominally have four coordinates even - # though the space is of dimension three only. The vector space V - # is supposed to hold the entire long vector, and the subspace W - # of V will be spanned by the vectors that arise from symmetric - # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3. - field = basis[0].base_ring() - dimension = basis[0].nrows() - - V = VectorSpace(field, dimension**2) - W = V.span_of_basis( _mat2vec(s) for s in basis ) - n = len(basis) - mult_table = [[W.zero() for j in xrange(n)] for i in xrange(n)] - for i in xrange(n): - for j in xrange(n): - mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2 - mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry)) + def jordan_product(X,Y): + return (X*Y + Y*X)/2 - return mult_table + @staticmethod + def trace_inner_product(X,Y): + r""" + A trace inner-product for matrices that aren't embedded in the + reals. It takes MATRICES as arguments, not EJA elements. + SETUP:: - @staticmethod - def real_embed(M): - """ - Embed the matrix ``M`` into a space of real matrices. + sage: from mjo.eja.eja_algebra import (RealSymmetricEJA, + ....: ComplexHermitianEJA, + ....: QuaternionHermitianEJA, + ....: OctonionHermitianEJA) - The matrix ``M`` can have entries in any field at the moment: - the real numbers, complex numbers, or quaternions. And although - they are not a field, we can probably support octonions at some - point, too. This function returns a real matrix that "acts like" - the original with respect to matrix multiplication; i.e. + EXAMPLES:: - real_embed(M*N) = real_embed(M)*real_embed(N) + sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 - """ - raise NotImplementedError + :: + sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 - @staticmethod - def real_unembed(M): - """ - The inverse of :meth:`real_embed`. - """ - raise NotImplementedError + :: + sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 - @classmethod - def natural_inner_product(cls,X,Y): - Xu = cls.real_unembed(X) - Yu = cls.real_unembed(Y) - tr = (Xu*Yu).trace() + :: - if tr in RLF: - # It's real already. - return tr + sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 - # Otherwise, try the thing that works for complex numbers; and - # if that doesn't work, the thing that works for quaternions. - try: - return tr.vector()[0] # real part, imag part is index 1 - except AttributeError: - # A quaternions doesn't have a vector() method, but does - # have coefficient_tuple() method that returns the - # coefficients of 1, i, j, and k -- in that order. + """ + tr = (X*Y).trace() + if hasattr(tr, 'coefficient'): + # Works for octonions, and has to come first because they + # also have a "real()" method that doesn't return an + # element of the scalar ring. + return tr.coefficient(0) + elif hasattr(tr, 'coefficient_tuple'): + # Works for quaternions. return tr.coefficient_tuple()[0] + # Works for real and complex numbers. + return tr.real() -class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): - @staticmethod - def real_embed(M): - """ - The identity function, for embedding real matrices into real - matrices. - """ - return M - @staticmethod - def real_unembed(M): - """ - The identity function, for unembedding real matrices from real - matrices. - """ - return M + def __init__(self, matrix_space, **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().__init__(self._denormalized_basis(matrix_space), + self.jordan_product, + self.trace_inner_product, + field=matrix_space.base_ring(), + matrix_space=matrix_space, + **kwargs) -class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): + self.rank.set_cache(matrix_space.nrows()) + self.one.set_cache( self(matrix_space.one()) ) + +class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -1069,19 +1916,19 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): 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, AA) - Euclidean Jordan algebra of dimension 3 over Algebraic Real Field - sage: RealSymmetricEJA(2, RR) + sage: RealSymmetricEJA(2, field=RDF, check_axioms=True) + Euclidean Jordan algebra of dimension 3 over Real Double Field + sage: RealSymmetricEJA(2, field=RR, check_axioms=True) Euclidean Jordan algebra of dimension 3 over Real Field with 53 bits of precision @@ -1090,7 +1937,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): The dimension of this algebra is `(n^2 + n) / 2`:: sage: set_random_seed() - sage: n_max = RealSymmetricEJA._max_test_case_size() + sage: n_max = RealSymmetricEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) sage: J = RealSymmetricEJA(n) sage: J.dimension() == (n^2 + n)/2 @@ -1101,9 +1948,9 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): sage: set_random_seed() sage: J = RealSymmetricEJA.random_instance() sage: x,y = J.random_elements(2) - sage: actual = (x*y).natural_representation() - sage: X = x.natural_representation() - sage: Y = y.natural_representation() + sage: actual = (x*y).to_matrix() + sage: X = x.to_matrix() + sage: Y = y.to_matrix() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True @@ -1115,262 +1962,168 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): sage: RealSymmetricEJA(3, prefix='q').gens() (q0, q1, q2, q3, q4, q5) - Our natural basis is normalized with respect to the natural inner - product unless we specify otherwise:: - - sage: set_random_seed() - sage: J = RealSymmetricEJA.random_instance() - sage: all( b.norm() == 1 for b in J.gens() ) - True - - Since our natural basis is normalized with respect to the natural - inner product, and since we know that this algebra is an EJA, any - left-multiplication operator's matrix will be symmetric because - natural->EJA basis representation is an isometry and within the EJA - the operator is self-adjoint by the Jordan axiom:: + We can construct the (trivial) algebra of rank zero:: - sage: set_random_seed() - sage: x = RealSymmetricEJA.random_instance().random_element() - sage: x.operator().matrix().is_symmetric() - True + sage: RealSymmetricEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ + @staticmethod + def _max_random_instance_size(): + return 4 # Dimension 10 + @classmethod - def _denormalized_basis(cls, n, field): + def random_instance(cls, **kwargs): """ - Return a basis for the space of real symmetric n-by-n matrices. - - SETUP:: - - sage: from mjo.eja.eja_algebra import RealSymmetricEJA + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, **kwargs) - TESTS:: + 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 - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: B = RealSymmetricEJA._denormalized_basis(n,QQ) - sage: all( M.is_symmetric() for M in B) - True + A = MatrixSpace(field, n) + super().__init__(A, **kwargs) - """ - # The basis of symmetric matrices, as matrices, in their R^(n-by-n) - # coordinates. - S = [] - for i in xrange(n): - for j in xrange(i+1): - Eij = matrix(field, n, lambda k,l: k==i and l==j) - if i == j: - Sij = Eij - else: - Sij = Eij + Eij.transpose() - S.append(Sij) - return S - @staticmethod - def _max_test_case_size(): - return 4 # Dimension 10 - +class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): + """ + The rank-n simple EJA consisting of complex Hermitian n-by-n + matrices over the real numbers, the usual symmetric Jordan product, + and the real-part-of-trace inner product. It has dimension `n^2` over + the reals. - def __init__(self, n, field=QQ, **kwargs): - basis = self._denormalized_basis(n, field) - super(RealSymmetricEJA, self).__init__(field, basis, n, **kwargs) + SETUP:: + sage: from mjo.eja.eja_algebra import ComplexHermitianEJA -class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): - @staticmethod - def real_embed(M): - """ - Embed the n-by-n complex matrix ``M`` into the space of real - matrices of size 2n-by-2n via the map the sends each entry `z = a + - bi` to the block matrix ``[[a,b],[-b,a]]``. + EXAMPLES: - SETUP:: + In theory, our "field" can be any subfield of the reals, but we + can't use inexact real fields at the moment because SageMath + doesn't know how to convert their elements into complex numbers, + or even into algebraic reals:: + + sage: QQbar(RDF(1)) + Traceback (most recent call last): + ... + TypeError: Illegal initializer for algebraic number + sage: AA(RR(1)) + Traceback (most recent call last): + ... + TypeError: Illegal initializer for algebraic number + + This causes the following error when we try to scale a matrix of + complex numbers by an inexact real number:: + + sage: ComplexHermitianEJA(2,field=RR) + Traceback (most recent call last): + ... + TypeError: Unable to coerce entries (=(1.00000000000000, + -0.000000000000000)) to coefficients in Algebraic Real Field - sage: from mjo.eja.eja_algebra import \ - ....: ComplexMatrixEuclideanJordanAlgebra + TESTS: - EXAMPLES:: + The dimension of this algebra is `n^2`:: - sage: F = QuadraticField(-1, 'i') - sage: x1 = F(4 - 2*i) - sage: x2 = F(1 + 2*i) - sage: x3 = F(-i) - sage: x4 = F(6) - sage: M = matrix(F,2,[[x1,x2],[x3,x4]]) - sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M) - [ 4 -2| 1 2] - [ 2 4|-2 1] - [-----+-----] - [ 0 -1| 6 0] - [ 1 0| 0 6] + sage: set_random_seed() + sage: n_max = ComplexHermitianEJA._max_random_instance_size() + sage: n = ZZ.random_element(1, n_max) + sage: J = ComplexHermitianEJA(n) + sage: J.dimension() == n^2 + True - TESTS: + The Jordan multiplication is what we think it is:: - Embedding is a homomorphism (isomorphism, in fact):: + sage: set_random_seed() + sage: J = ComplexHermitianEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: actual = (x*y).to_matrix() + sage: X = x.to_matrix() + sage: Y = y.to_matrix() + sage: expected = (X*Y + Y*X)/2 + sage: actual == expected + True + sage: J(expected) == x*y + True - sage: set_random_seed() - sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size() - sage: n = ZZ.random_element(n_max) - sage: F = QuadraticField(-1, 'i') - sage: X = random_matrix(F, n) - sage: Y = random_matrix(F, n) - sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X) - sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y) - sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y) - sage: Xe*Ye == XYe - True + We can change the generator prefix:: - """ - n = M.nrows() - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") + sage: ComplexHermitianEJA(2, prefix='z').gens() + (z0, z1, z2, z3) - # We don't need any adjoined elements... - field = M.base_ring().base_ring() + We can construct the (trivial) algebra of rank zero:: - blocks = [] - for z in M.list(): - a = z.list()[0] # real part, I guess - b = z.list()[1] # imag part, I guess - blocks.append(matrix(field, 2, [[a,b],[-b,a]])) + sage: ComplexHermitianEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ + 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 - return matrix.block(field, n, blocks) + from mjo.hurwitz import ComplexMatrixAlgebra + A = ComplexMatrixAlgebra(n, scalars=field) + super().__init__(A, **kwargs) @staticmethod - def real_unembed(M): - """ - The inverse of _embed_complex_matrix(). - - SETUP:: - - sage: from mjo.eja.eja_algebra import \ - ....: ComplexMatrixEuclideanJordanAlgebra - - EXAMPLES:: - - sage: A = matrix(QQ,[ [ 1, 2, 3, 4], - ....: [-2, 1, -4, 3], - ....: [ 9, 10, 11, 12], - ....: [-10, 9, -12, 11] ]) - sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A) - [ 2*i + 1 4*i + 3] - [ 10*i + 9 12*i + 11] - - TESTS: - - Unembedding is the inverse of embedding:: - - sage: set_random_seed() - sage: F = QuadraticField(-1, 'i') - sage: M = random_matrix(F, 3) - sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M) - sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M - True - - """ - n = ZZ(M.nrows()) - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - if not n.mod(2).is_zero(): - raise ValueError("the matrix 'M' must be a complex embedding") - - # If "M" was normalized, its base ring might have roots - # adjoined and they can stick around after unembedding. - field = M.base_ring() - R = PolynomialRing(field, 'z') - z = R.gen() - F = field.extension(z**2 + 1, 'i', embedding=CLF(-1).sqrt()) - i = F.gen() - - # Go top-left to bottom-right (reading order), converting every - # 2-by-2 block we see to a single complex element. - elements = [] - for k in xrange(n/2): - for j in xrange(n/2): - submat = M[2*k:2*k+2,2*j:2*j+2] - if submat[0,0] != submat[1,1]: - raise ValueError('bad on-diagonal submatrix') - if submat[0,1] != -submat[1,0]: - raise ValueError('bad off-diagonal submatrix') - z = submat[0,0] + submat[0,1]*i - elements.append(z) - - return matrix(F, n/2, elements) - + def _max_random_instance_size(): + return 3 # Dimension 9 @classmethod - def natural_inner_product(cls,X,Y): + def random_instance(cls, **kwargs): """ - Compute a natural inner product in this algebra directly from - its real embedding. - - SETUP:: - - sage: from mjo.eja.eja_algebra import ComplexHermitianEJA - - TESTS: - - This gives the same answer as the slow, default method implemented - in :class:`MatrixEuclideanJordanAlgebra`:: - - sage: set_random_seed() - sage: J = ComplexHermitianEJA.random_instance() - sage: x,y = J.random_elements(2) - sage: Xe = x.natural_representation() - sage: Ye = y.natural_representation() - sage: X = ComplexHermitianEJA.real_unembed(Xe) - sage: Y = ComplexHermitianEJA.real_unembed(Ye) - sage: expected = (X*Y).trace().vector()[0] - sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye) - sage: actual == expected - True - + Return a random instance of this type of algebra. """ - return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2 + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, **kwargs) -class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): - """ - The rank-n simple EJA consisting of complex Hermitian n-by-n - matrices over the real numbers, the usual symmetric Jordan product, - and the real-part-of-trace inner product. It has dimension `n^2` over +class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): + r""" + The rank-n simple EJA consisting of self-adjoint n-by-n quaternion + matrices, the usual symmetric Jordan product, and the + real-part-of-trace inner product. It has dimension `2n^2 - n` over the reals. SETUP:: - sage: from mjo.eja.eja_algebra import ComplexHermitianEJA + sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA EXAMPLES: In theory, our "field" can be any subfield of the reals:: - sage: ComplexHermitianEJA(2, AA) - Euclidean Jordan algebra of dimension 4 over Algebraic Real Field - sage: ComplexHermitianEJA(2, RR) - Euclidean Jordan algebra of dimension 4 over Real Field with + sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True) + Euclidean Jordan algebra of dimension 6 over Real Double Field + sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True) + Euclidean Jordan algebra of dimension 6 over Real Field with 53 bits of precision TESTS: - The dimension of this algebra is `n^2`:: + The dimension of this algebra is `2*n^2 - n`:: sage: set_random_seed() - sage: n_max = ComplexHermitianEJA._max_test_case_size() + sage: n_max = QuaternionHermitianEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) - sage: J = ComplexHermitianEJA(n) - sage: J.dimension() == n^2 + sage: J = QuaternionHermitianEJA(n) + sage: J.dimension() == 2*(n^2) - n True The Jordan multiplication is what we think it is:: sage: set_random_seed() - sage: J = ComplexHermitianEJA.random_instance() + sage: J = QuaternionHermitianEJA.random_instance() sage: x,y = J.random_elements(2) - sage: actual = (x*y).natural_representation() - sage: X = x.natural_representation() - sage: Y = y.natural_representation() + sage: actual = (x*y).to_matrix() + sage: X = x.to_matrix() + sage: Y = y.to_matrix() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True @@ -1379,390 +2132,439 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): We can change the generator prefix:: - sage: ComplexHermitianEJA(2, prefix='z').gens() - (z0, z1, z2, z3) + sage: QuaternionHermitianEJA(2, prefix='a').gens() + (a0, a1, a2, a3, a4, a5) - Our natural basis is normalized with respect to the natural inner - product unless we specify otherwise:: + We can construct the (trivial) algebra of rank zero:: - sage: set_random_seed() - sage: J = ComplexHermitianEJA.random_instance() - sage: all( b.norm() == 1 for b in J.gens() ) - True + sage: QuaternionHermitianEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field - Since our natural basis is normalized with respect to the natural - inner product, and since we know that this algebra is an EJA, any - left-multiplication operator's matrix will be symmetric because - natural->EJA basis representation is an isometry and within the EJA - the operator is self-adjoint by the Jordan axiom:: + """ + 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 - sage: set_random_seed() - sage: x = ComplexHermitianEJA.random_instance().random_element() - sage: x.operator().matrix().is_symmetric() - True + from mjo.hurwitz import QuaternionMatrixAlgebra + A = QuaternionMatrixAlgebra(n, scalars=field) + super().__init__(A, **kwargs) - """ - @classmethod - def _denormalized_basis(cls, n, field): + @staticmethod + def _max_random_instance_size(): + r""" + The maximum rank of a random QuaternionHermitianEJA. """ - Returns a basis for the space of complex Hermitian n-by-n matrices. + return 2 # Dimension 6 - Why do we embed these? Basically, because all of numerical linear - algebra assumes that you're working with vectors consisting of `n` - entries from a field and scalars from the same field. There's no way - to tell SageMath that (for example) the vectors contain complex - numbers, while the scalar field is real. + @classmethod + def random_instance(cls, **kwargs): + """ + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, **kwargs) - SETUP:: +class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): + r""" + SETUP:: - sage: from mjo.eja.eja_algebra import ComplexHermitianEJA + sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA, + ....: OctonionHermitianEJA) + sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra - TESTS:: + EXAMPLES: - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: field = QuadraticField(2, 'sqrt2') - sage: B = ComplexHermitianEJA._denormalized_basis(n, field) - sage: all( M.is_symmetric() for M in B) - True + The 3-by-3 algebra satisfies the axioms of an EJA:: + + sage: OctonionHermitianEJA(3, # long time + ....: field=QQ, # long time + ....: orthonormalize=False, # long time + ....: check_axioms=True) # long time + Euclidean Jordan algebra of dimension 27 over Rational Field + + After a change-of-basis, the 2-by-2 algebra has the same + multiplication table as the ten-dimensional Jordan spin algebra:: + + sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ) + sage: b = OctonionHermitianEJA._denormalized_basis(A) + sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],) + sage: jp = OctonionHermitianEJA.jordan_product + sage: ip = OctonionHermitianEJA.trace_inner_product + sage: J = FiniteDimensionalEJA(basis, + ....: jp, + ....: ip, + ....: field=QQ, + ....: orthonormalize=False) + sage: J.multiplication_table() + +----++----+----+----+----+----+----+----+----+----+----+ + | * || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 | + +====++====+====+====+====+====+====+====+====+====+====+ + | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b1 || b1 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b2 || b2 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b3 || b3 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b4 || b4 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b5 || b5 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b6 || b6 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b7 || b7 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b8 || b8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b9 || b9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | + +----++----+----+----+----+----+----+----+----+----+----+ - """ - R = PolynomialRing(field, 'z') - z = R.gen() - F = field.extension(z**2 + 1, 'I') - I = F.gen() + TESTS: - # This is like the symmetric case, but we need to be careful: - # - # * We want conjugate-symmetry, not just symmetry. - # * The diagonal will (as a result) be real. - # - S = [] - for i in xrange(n): - for j in xrange(i+1): - Eij = matrix(F, n, lambda k,l: k==i and l==j) - if i == j: - Sij = cls.real_embed(Eij) - S.append(Sij) - else: - # The second one has a minus because it's conjugated. - Sij_real = cls.real_embed(Eij + Eij.transpose()) - S.append(Sij_real) - Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose()) - S.append(Sij_imag) + We can actually construct the 27-dimensional Albert algebra, + and we get the right unit element if we recompute it:: + + sage: J = OctonionHermitianEJA(3, # long time + ....: field=QQ, # long time + ....: orthonormalize=False) # long time + sage: J.one.clear_cache() # long time + sage: J.one() # long time + b0 + b9 + b26 + sage: J.one().to_matrix() # long time + +----+----+----+ + | e0 | 0 | 0 | + +----+----+----+ + | 0 | e0 | 0 | + +----+----+----+ + | 0 | 0 | e0 | + +----+----+----+ + + The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan + spin algebra, but just to be sure, we recompute its rank:: + + sage: J = OctonionHermitianEJA(2, # long time + ....: field=QQ, # long time + ....: orthonormalize=False) # long time + sage: J.rank.clear_cache() # long time + sage: J.rank() # long time + 2 - # Since we embedded these, we can drop back to the "field" that we - # started with instead of the complex extension "F". - return ( s.change_ring(field) for s in S ) + """ + @staticmethod + def _max_random_instance_size(): + r""" + The maximum rank of a random QuaternionHermitianEJA. + """ + return 1 # Dimension 1 + @classmethod + def random_instance(cls, **kwargs): + """ + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, **kwargs) - def __init__(self, n, field=QQ, **kwargs): - basis = self._denormalized_basis(n,field) - super(ComplexHermitianEJA,self).__init__(field, basis, n, **kwargs) + def __init__(self, n, field=AA, **kwargs): + if n > 3: + # Otherwise we don't get an EJA. + raise ValueError("n cannot exceed 3") + # 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 -class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): - @staticmethod - def real_embed(M): - """ - Embed the n-by-n quaternion matrix ``M`` into the space of real - matrices of size 4n-by-4n by first sending each quaternion entry `z - = a + bi + cj + dk` to the block-complex matrix ``[[a + bi, - c+di],[-c + di, a-bi]]`, and then embedding those into a real - matrix. + from mjo.hurwitz import OctonionMatrixAlgebra + A = OctonionMatrixAlgebra(n, scalars=field) + super().__init__(A, **kwargs) - SETUP:: - sage: from mjo.eja.eja_algebra import \ - ....: QuaternionMatrixEuclideanJordanAlgebra +class AlbertEJA(OctonionHermitianEJA): + r""" + The Albert algebra is the algebra of three-by-three Hermitian + matrices whose entries are octonions. - EXAMPLES:: + SETUP:: - sage: Q = QuaternionAlgebra(QQ,-1,-1) - sage: i,j,k = Q.gens() - sage: x = 1 + 2*i + 3*j + 4*k - sage: M = matrix(Q, 1, [[x]]) - sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M) - [ 1 2 3 4] - [-2 1 -4 3] - [-3 4 1 -2] - [-4 -3 2 1] + sage: from mjo.eja.eja_algebra import AlbertEJA - Embedding is a homomorphism (isomorphism, in fact):: + EXAMPLES:: - sage: set_random_seed() - sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size() - sage: n = ZZ.random_element(n_max) - sage: Q = QuaternionAlgebra(QQ,-1,-1) - sage: X = random_matrix(Q, n) - sage: Y = random_matrix(Q, n) - sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X) - sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y) - sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y) - sage: Xe*Ye == XYe - True + sage: AlbertEJA(field=QQ, orthonormalize=False) + Euclidean Jordan algebra of dimension 27 over Rational Field + sage: AlbertEJA() # long time + Euclidean Jordan algebra of dimension 27 over Algebraic Real Field - """ - quaternions = M.base_ring() - n = M.nrows() - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") + """ + def __init__(self, *args, **kwargs): + super().__init__(3, *args, **kwargs) - F = QuadraticField(-1, 'i') - i = F.gen() - blocks = [] - for z in M.list(): - t = z.coefficient_tuple() - a = t[0] - b = t[1] - c = t[2] - d = t[3] - cplxM = matrix(F, 2, [[ a + b*i, c + d*i], - [-c + d*i, a - b*i]]) - realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM) - blocks.append(realM) +class HadamardEJA(RationalBasisEJA, ConcreteEJA): + """ + Return the Euclidean Jordan algebra on `R^n` with the Hadamard + (pointwise real-number multiplication) Jordan product and the + usual inner-product. - # We should have real entries by now, so use the realest field - # we've got for the return value. - return matrix.block(quaternions.base_ring(), n, blocks) + This is nothing more than the Cartesian product of ``n`` copies of + the one-dimensional Jordan spin algebra, and is the most common + example of a non-simple Euclidean Jordan algebra. + SETUP:: + sage: from mjo.eja.eja_algebra import HadamardEJA - @staticmethod - def real_unembed(M): - """ - The inverse of _embed_quaternion_matrix(). - - SETUP:: + EXAMPLES: - sage: from mjo.eja.eja_algebra import \ - ....: QuaternionMatrixEuclideanJordanAlgebra + This multiplication table can be verified by hand:: - EXAMPLES:: + sage: J = HadamardEJA(3) + sage: b0,b1,b2 = J.gens() + sage: b0*b0 + b0 + sage: b0*b1 + 0 + sage: b0*b2 + 0 + sage: b1*b1 + b1 + sage: b1*b2 + 0 + sage: b2*b2 + b2 - sage: M = matrix(QQ, [[ 1, 2, 3, 4], - ....: [-2, 1, -4, 3], - ....: [-3, 4, 1, -2], - ....: [-4, -3, 2, 1]]) - sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M) - [1 + 2*i + 3*j + 4*k] + TESTS: - TESTS: + We can change the generator prefix:: - Unembedding is the inverse of embedding:: + sage: HadamardEJA(3, prefix='r').gens() + (r0, r1, r2) + """ + def __init__(self, n, field=AA, **kwargs): + MS = MatrixSpace(field, n, 1) - sage: set_random_seed() - sage: Q = QuaternionAlgebra(QQ, -1, -1) - sage: M = random_matrix(Q, 3) - sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M) - sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M - True + if n == 0: + jordan_product = lambda x,y: x + inner_product = lambda x,y: x + else: + def jordan_product(x,y): + return MS( xi*yi for (xi,yi) in zip(x,y) ) + + def inner_product(x,y): + return (x.T*y)[0,0] + + # New defaults for keyword arguments. Don't orthonormalize + # because our basis is already orthonormal with respect to our + # inner-product. Don't check the axioms, because we know this + # is a valid EJA... but do double-check if the user passes + # check_axioms=True. Note: we DON'T override the "check_field" + # default here, because the user can pass in a field! + if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False + if "check_axioms" not in kwargs: kwargs["check_axioms"] = False + + column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() ) + super().__init__(column_basis, + jordan_product, + inner_product, + field=field, + matrix_space=MS, + associative=True, + **kwargs) + self.rank.set_cache(n) + + self.one.set_cache( self.sum(self.gens()) ) + @staticmethod + def _max_random_instance_size(): + r""" + The maximum dimension of a random HadamardEJA. """ - n = ZZ(M.nrows()) - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - if not n.mod(4).is_zero(): - raise ValueError("the matrix 'M' must be a quaternion embedding") - - # Use the base ring of the matrix to ensure that its entries can be - # multiplied by elements of the quaternion algebra. - field = M.base_ring() - Q = QuaternionAlgebra(field,-1,-1) - i,j,k = Q.gens() - - # Go top-left to bottom-right (reading order), converting every - # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1 - # quaternion block. - elements = [] - for l in xrange(n/4): - for m in xrange(n/4): - submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed( - M[4*l:4*l+4,4*m:4*m+4] ) - if submat[0,0] != submat[1,1].conjugate(): - raise ValueError('bad on-diagonal submatrix') - if submat[0,1] != -submat[1,0].conjugate(): - raise ValueError('bad off-diagonal submatrix') - z = submat[0,0].vector()[0] # real part - z += submat[0,0].vector()[1]*i # imag part - z += submat[0,1].vector()[0]*j # real part - z += submat[0,1].vector()[1]*k # imag part - elements.append(z) - - return matrix(Q, n/4, elements) - + return 5 @classmethod - def natural_inner_product(cls,X,Y): + def random_instance(cls, **kwargs): """ - Compute a natural inner product in this algebra directly from - its real embedding. - - SETUP:: - - sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA - - TESTS: - - This gives the same answer as the slow, default method implemented - in :class:`MatrixEuclideanJordanAlgebra`:: - - sage: set_random_seed() - sage: J = QuaternionHermitianEJA.random_instance() - sage: x,y = J.random_elements(2) - sage: Xe = x.natural_representation() - sage: Ye = y.natural_representation() - sage: X = QuaternionHermitianEJA.real_unembed(Xe) - sage: Y = QuaternionHermitianEJA.real_unembed(Ye) - sage: expected = (X*Y).trace().coefficient_tuple()[0] - sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye) - sage: actual == expected - True - + Return a random instance of this type of algebra. """ - return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4 + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, **kwargs) -class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, - KnownRankEJA): - """ - The rank-n simple EJA consisting of self-adjoint n-by-n quaternion - matrices, the usual symmetric Jordan product, and the - real-part-of-trace inner product. It has dimension `2n^2 - n` over - the reals. +class BilinearFormEJA(RationalBasisEJA, ConcreteEJA): + r""" + The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)`` + with the half-trace inner product and jordan product ``x*y = + (,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is + a symmetric positive-definite "bilinear form" matrix. Its + dimension is the size of `B`, and it has rank two in dimensions + larger than two. It reduces to the ``JordanSpinEJA`` when `B` is + the identity matrix of order ``n``. + + We insist that the one-by-one upper-left identity block of `B` be + passed in as well so that we can be passed a matrix of size zero + to construct a trivial algebra. SETUP:: - sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA + sage: from mjo.eja.eja_algebra import (BilinearFormEJA, + ....: JordanSpinEJA) EXAMPLES: - In theory, our "field" can be any subfield of the reals:: - - sage: QuaternionHermitianEJA(2, AA) - Euclidean Jordan algebra of dimension 6 over Algebraic Real Field - sage: QuaternionHermitianEJA(2, RR) - Euclidean Jordan algebra of dimension 6 over Real Field with - 53 bits of precision - - TESTS: + When no bilinear form is specified, the identity matrix is used, + and the resulting algebra is the Jordan spin algebra:: - The dimension of this algebra is `2*n^2 - n`:: - - sage: set_random_seed() - sage: n_max = QuaternionHermitianEJA._max_test_case_size() - sage: n = ZZ.random_element(1, n_max) - sage: J = QuaternionHermitianEJA(n) - sage: J.dimension() == 2*(n^2) - n + sage: B = matrix.identity(AA,3) + sage: J0 = BilinearFormEJA(B) + sage: J1 = JordanSpinEJA(3) + sage: J0.multiplication_table() == J0.multiplication_table() True - The Jordan multiplication is what we think it is:: - - sage: set_random_seed() - sage: J = QuaternionHermitianEJA.random_instance() - sage: x,y = J.random_elements(2) - sage: actual = (x*y).natural_representation() - sage: X = x.natural_representation() - sage: Y = y.natural_representation() - sage: expected = (X*Y + Y*X)/2 - sage: actual == expected - True - sage: J(expected) == x*y - True + An error is raised if the matrix `B` does not correspond to a + positive-definite bilinear form:: - We can change the generator prefix:: + sage: B = matrix.random(QQ,2,3) + sage: J = BilinearFormEJA(B) + Traceback (most recent call last): + ... + ValueError: bilinear form is not positive-definite + sage: B = matrix.zero(QQ,3) + sage: J = BilinearFormEJA(B) + Traceback (most recent call last): + ... + ValueError: bilinear form is not positive-definite - sage: QuaternionHermitianEJA(2, prefix='a').gens() - (a0, a1, a2, a3, a4, a5) + TESTS: - Our natural basis is normalized with respect to the natural inner - product unless we specify otherwise:: + We can create a zero-dimensional algebra:: - sage: set_random_seed() - sage: J = QuaternionHermitianEJA.random_instance() - sage: all( b.norm() == 1 for b in J.gens() ) - True + sage: B = matrix.identity(AA,0) + sage: J = BilinearFormEJA(B) + sage: J.basis() + Finite family {} - Since our natural basis is normalized with respect to the natural - inner product, and since we know that this algebra is an EJA, any - left-multiplication operator's matrix will be symmetric because - natural->EJA basis representation is an isometry and within the EJA - the operator is self-adjoint by the Jordan axiom:: + We can check the multiplication condition given in the Jordan, von + Neumann, and Wigner paper (and also discussed on my "On the + symmetry..." paper). Note that this relies heavily on the standard + choice of basis, as does anything utilizing the bilinear form + matrix. We opt not to orthonormalize the basis, because if we + did, we would have to normalize the `s_{i}` in a similar manner:: sage: set_random_seed() - sage: x = QuaternionHermitianEJA.random_instance().random_element() - sage: x.operator().matrix().is_symmetric() + sage: n = ZZ.random_element(5) + sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular') + sage: B11 = matrix.identity(QQ,1) + sage: B22 = M.transpose()*M + sage: B = block_matrix(2,2,[ [B11,0 ], + ....: [0, B22 ] ]) + sage: J = BilinearFormEJA(B, orthonormalize=False) + sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis() + sage: V = J.vector_space() + sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() ) + ....: for ei in eis ] + sage: actual = [ sis[i]*sis[j] + ....: for i in range(n-1) + ....: for j in range(n-1) ] + sage: expected = [ J.one() if i == j else J.zero() + ....: for i in range(n-1) + ....: for j in range(n-1) ] + sage: actual == expected True """ - @classmethod - def _denormalized_basis(cls, n, field): - """ - Returns a basis for the space of quaternion Hermitian n-by-n matrices. - - Why do we embed these? Basically, because all of numerical - linear algebra assumes that you're working with vectors consisting - of `n` entries from a field and scalars from the same field. There's - no way to tell SageMath that (for example) the vectors contain - complex numbers, while the scalar field is real. - - SETUP:: - - sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA - - TESTS:: - - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ) - sage: all( M.is_symmetric() for M in B ) - True + 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... + if ("check_axioms" not in kwargs) or kwargs["check_axioms"]: + if not B.is_positive_definite(): + raise ValueError("bilinear form is not positive-definite") + + # However, all of the other data for this EJA is computed + # by us in manner that guarantees the axioms are + # satisfied. So, again, unless we are specifically asked to + # verify things, we'll skip the rest of the checks. + if "check_axioms" not in kwargs: kwargs["check_axioms"] = False + + n = B.nrows() + MS = MatrixSpace(field, n, 1) + + def inner_product(x,y): + return (y.T*B*x)[0,0] + + def jordan_product(x,y): + x0 = x[0,0] + xbar = x[1:,0] + y0 = y[0,0] + ybar = y[1:,0] + z0 = inner_product(y,x) + zbar = y0*xbar + x0*ybar + return MS([z0] + zbar.list()) + + column_basis = tuple( MS(b) 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, + matrix_space=MS, + associative=associative, + **kwargs) + + # The rank of this algebra is two, unless we're in a + # one-dimensional ambient space (because the rank is bounded + # by the ambient dimension). + self.rank.set_cache(min(n,2)) + if n == 0: + self.one.set_cache( self.zero() ) + else: + self.one.set_cache( self.monomial(0) ) + @staticmethod + def _max_random_instance_size(): + r""" + The maximum dimension of a random BilinearFormEJA. """ - Q = QuaternionAlgebra(QQ,-1,-1) - I,J,K = Q.gens() + return 5 - # This is like the symmetric case, but we need to be careful: - # - # * We want conjugate-symmetry, not just symmetry. - # * The diagonal will (as a result) be real. - # - S = [] - for i in xrange(n): - for j in xrange(i+1): - Eij = matrix(Q, n, lambda k,l: k==i and l==j) - if i == j: - Sij = cls.real_embed(Eij) - S.append(Sij) - else: - # The second, third, and fourth ones have a minus - # because they're conjugated. - Sij_real = cls.real_embed(Eij + Eij.transpose()) - S.append(Sij_real) - Sij_I = cls.real_embed(I*Eij - I*Eij.transpose()) - S.append(Sij_I) - Sij_J = cls.real_embed(J*Eij - J*Eij.transpose()) - S.append(Sij_J) - Sij_K = cls.real_embed(K*Eij - K*Eij.transpose()) - S.append(Sij_K) + @classmethod + def random_instance(cls, **kwargs): + """ + Return a random instance of this algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + if n.is_zero(): + B = matrix.identity(ZZ, n) + return cls(B, **kwargs) - # Since we embedded these, we can drop back to the "field" that we - # started with instead of the quaternion algebra "Q". - return ( s.change_ring(field) for s in S ) + B11 = matrix.identity(ZZ, 1) + M = matrix.random(ZZ, n-1) + I = matrix.identity(ZZ, n-1) + alpha = ZZ.zero() + while alpha.is_zero(): + alpha = ZZ.random_element().abs() + B22 = M.transpose()*M + alpha*I + from sage.matrix.special import block_matrix + B = block_matrix(2,2, [ [B11, ZZ(0) ], + [ZZ(0), B22 ] ]) - def __init__(self, n, field=QQ, **kwargs): - basis = self._denormalized_basis(n,field) - super(QuaternionHermitianEJA,self).__init__(field, basis, n, **kwargs) + return cls(B, **kwargs) -class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): +class JordanSpinEJA(BilinearFormEJA): """ The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)`` with the usual inner product and jordan product ``x*y = - (, x0*y_bar + y0*x_bar)``. It has dimension `n` over + (, x0*y_bar + y0*x_bar)``. It has dimension `n` over the reals. SETUP:: @@ -1774,20 +2576,20 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): 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:: @@ -1795,56 +2597,51 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): sage: JordanSpinEJA(2, prefix='B').gens() (B0, B1) - """ - def __init__(self, n, field=QQ, **kwargs): - V = VectorSpace(field, n) - mult_table = [[V.zero() for j in xrange(n)] for i in xrange(n)] - for i in xrange(n): - for j in xrange(n): - x = V.gen(i) - y = V.gen(j) - x0 = x[0] - xbar = x[1:] - y0 = y[0] - ybar = y[1:] - # z = x*y - z0 = x.inner_product(y) - zbar = y0*xbar + x0*ybar - z = V([z0] + zbar.list()) - mult_table[i][j] = z - - # The rank of the spin algebra is two, unless we're in a - # one-dimensional ambient space (because the rank is bounded by - # the ambient dimension). - fdeja = super(JordanSpinEJA, self) - return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs) - - def inner_product(self, x, y): - """ - Faster to reimplement than to use natural representations. - - SETUP:: - - sage: from mjo.eja.eja_algebra import JordanSpinEJA - - TESTS: + TESTS: - Ensure that this is the usual inner product for the algebras - over `R^n`:: + Ensure that we have the usual inner product on `R^n`:: sage: set_random_seed() sage: J = JordanSpinEJA.random_instance() sage: x,y = J.random_elements(2) - sage: X = x.natural_representation() - sage: Y = y.natural_representation() - sage: x.inner_product(y) == J.natural_inner_product(X,Y) + sage: actual = x.inner_product(y) + sage: expected = x.to_vector().inner_product(y.to_vector()) + sage: actual == expected True + """ + 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) + + # Don't orthonormalize because our basis is already + # orthonormal with respect to our inner-product. + if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False + + # But also don't pass check_field=False here, because the user + # can pass in a field! + super().__init__(B, *args, **kwargs) + + @staticmethod + def _max_random_instance_size(): + r""" + The maximum dimension of a random JordanSpinEJA. """ - return x.to_vector().inner_product(y.to_vector()) + return 5 + @classmethod + def random_instance(cls, **kwargs): + """ + Return a random instance of this type of algebra. -class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): + Needed here to override the implementation for ``BilinearFormEJA``. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, **kwargs) + + +class TrivialEJA(RationalBasisEJA, ConcreteEJA): """ The trivial Euclidean Jordan algebra consisting of only a zero element. @@ -1868,14 +2665,591 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): sage: J.one().norm() 0 sage: J.one().subalgebra_generated_by() - Euclidean Jordan algebra of dimension 0 over Rational Field + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field sage: J.rank() 0 """ - def __init__(self, field=QQ, **kwargs): - mult_table = [] - fdeja = super(TrivialEJA, self) + def __init__(self, field=AA, **kwargs): + jordan_product = lambda x,y: x + inner_product = lambda x,y: field.zero() + basis = () + MS = MatrixSpace(field,0) + + # New defaults for keyword arguments + if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False + if "check_axioms" not in kwargs: kwargs["check_axioms"] = False + + super().__init__(basis, + jordan_product, + inner_product, + associative=True, + field=field, + matrix_space=MS, + **kwargs) + # The rank is zero using my definition, namely the dimension of the # largest subalgebra generated by any element. - return fdeja.__init__(field, mult_table, rank=0, **kwargs) + self.rank.set_cache(0) + self.one.set_cache( self.zero() ) + + @classmethod + def random_instance(cls, **kwargs): + # We don't take a "size" argument so the superclass method is + # inappropriate for us. + return cls(**kwargs) + + +class CartesianProductEJA(FiniteDimensionalEJA): + r""" + 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, + ....: JordanSpinEJA, + ....: RealSymmetricEJA) + + 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 = 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() + 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: + + All factors must share the same base field:: + + sage: J1 = HadamardEJA(2, field=QQ) + sage: J2 = RealSymmetricEJA(2) + sage: CartesianProductEJA((J1,J2)) + Traceback (most recent call last): + ... + 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 + + """ + Element = FiniteDimensionalEJAElement + + + def __init__(self, factors, **kwargs): + m = len(factors) + if m == 0: + return TrivialEJA() + + self._sets = factors + + field = factors[0].base_ring() + if not all( J.base_ring() == field for J in factors ): + raise ValueError("all factors must share the same base field") + + associative = all( f.is_associative() for f in factors ) + + # Compute my matrix space. This category isn't perfect, but + # is good enough for what we need to do. + MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis() + MS_cat = MS_cat.Unital().CartesianProducts() + MS_factors = tuple( J.matrix_space() for J in factors ) + from sage.sets.cartesian_product import CartesianProduct + MS = CartesianProduct(MS_factors, MS_cat) + + 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) + ) + + # 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, + matrix_space=MS, + orthonormalize=False, + associative=associative, + cartesian_product=True, + check_field=False, + check_axioms=False) + + self.rank.set_cache(sum(J.rank() for J in factors)) + ones = tuple(J.one().to_matrix() for J in factors) + self.one.set_cache(self(ones)) + + def cartesian_factors(self): + # Copy/pasted from CombinatorialFreeModule_CartesianProduct. + return self._sets + + def cartesian_factor(self, i): + r""" + Return the ``i``th factor of this algebra. + """ + return self._sets[i] + + def _repr_(self): + # Copy/pasted from CombinatorialFreeModule_CartesianProduct. + from sage.categories.cartesian_product import cartesian_product + return cartesian_product.symbol.join("%s" % factor + for factor in self._sets) + + def matrix_space(self): + r""" + Return the space that our matrix basis lives in as a Cartesian + product. + + We don't simply use the ``cartesian_product()`` functor here + because it acts differently on SageMath MatrixSpaces and our + custom MatrixAlgebras, which are CombinatorialFreeModules. We + always want the result to be represented (and indexed) as + an ordered tuple. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, + ....: HadamardEJA, + ....: OctonionHermitianEJA, + ....: RealSymmetricEJA) + + 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) + + :: + + sage: J1 = ComplexHermitianEJA(1) + sage: J2 = ComplexHermitianEJA(1) + sage: J = cartesian_product([J1,J2]) + sage: J.one().to_matrix()[0] + +---+ + | 1 | + +---+ + sage: J.one().to_matrix()[1] + +---+ + | 1 | + +---+ + + :: + + sage: J1 = OctonionHermitianEJA(1) + sage: J2 = OctonionHermitianEJA(1) + sage: J = cartesian_product([J1,J2]) + sage: J.one().to_matrix()[0] + +----+ + | e0 | + +----+ + sage: J.one().to_matrix()[1] + +----+ + | e0 | + +----+ + + """ + return super().matrix_space() + + + @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 (HadamardEJA, + ....: JordanSpinEJA, + ....: OctonionHermitianEJA, + ....: 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 + + TESTS: + + The ``cartesian_product()`` function only uses the first factor to + decide where the result will live; thus we have to be careful to + check that all factors do indeed have a `_rational_algebra` member + before we try to access it:: + + sage: J1 = OctonionHermitianEJA(1) # no rational basis + sage: J2 = HadamardEJA(2) + sage: cartesian_product([J1,J2]) + Euclidean Jordan algebra of dimension 1 over Algebraic Real Field + (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field + sage: cartesian_product([J2,J1]) + Euclidean Jordan algebra of dimension 2 over Algebraic Real Field + (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field + + """ + def __init__(self, algebras, **kwargs): + CartesianProductEJA.__init__(self, algebras, **kwargs) + + self._rational_algebra = None + if self.vector_space().base_field() is not QQ: + if all( hasattr(r, "_rational_algebra") for r in algebras ): + 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