"""
-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 `<x*y,z> = <y,x*z>` for all elements
+ `x,y,z` in the algebra.
+
+Every Euclidean Jordan algebra is formally-real: for any two elements
+`x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
+0`. Conversely, every finite-dimensional formally-real Jordan algebra
+can be made into a Euclidean Jordan algebra with an appropriate choice
+of inner-product.
+
+Formally-real Jordan algebras were originally studied as a framework
+for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
+symmetric cone optimization, since every symmetric cone arises as the
+cone of squares in some Euclidean Jordan algebra.
+
+It is known that every Euclidean Jordan algebra decomposes into an
+orthogonal direct sum (essentially, a Cartesian product) of simple
+algebras, and that moreover, up to Jordan-algebra isomorphism, there
+are only five families of simple algebras. We provide constructions
+for these simple algebras:
+
+ * :class:`BilinearFormEJA`
+ * :class:`RealSymmetricEJA`
+ * :class:`ComplexHermitianEJA`
+ * :class:`QuaternionHermitianEJA`
+
+Missing from this list is the algebra of three-by-three octononion
+Hermitian matrices, as there is (as of yet) no implementation of the
+octonions in SageMath. In addition to these, we provide two other
+example constructions,
+
+ * :class:`HadamardEJA`
+ * :class:`TrivialEJA`
+
+The Jordan spin algebra is a bilinear form algebra where the bilinear
+form is the identity. The Hadamard EJA is simply a Cartesian product
+of one-dimensional spin algebras. And last but not least, the trivial
+EJA is exactly what you think. Cartesian products of these are also
+supported using the usual ``cartesian_product()`` function; as a
+result, we support (up to isomorphism) all Euclidean Jordan algebras
+that don't involve octonions.
+
+SETUP::
+
+ sage: from mjo.eja.eja_algebra import random_eja
+
+EXAMPLES::
+
+ sage: random_eja()
+ Euclidean Jordan algebra of dimension...
"""
from itertools import repeat
from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
from sage.categories.magmatic_algebras import MagmaticAlgebras
+from sage.categories.sets_cat import cartesian_product
from sage.combinat.free_module import CombinatorialFreeModule
from sage.matrix.constructor import matrix
from sage.matrix.matrix_space import MatrixSpace
from sage.misc.cachefunc import cached_method
-from sage.misc.lazy_import import lazy_import
-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, AA, QQbar, RR, RLF, CLF,
PolynomialRing,
QuadraticField)
-from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
-lazy_import('mjo.eja.eja_subalgebra',
- 'FiniteDimensionalEuclideanJordanSubalgebra')
-from mjo.eja.eja_utils import _mat2vec
+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.
+
+ - ``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,
+ basis,
+ jordan_product,
+ inner_product,
+ field=AA,
+ 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 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...
+ self._matrix_basis = basis
+
+ # 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")
-class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
def _coerce_map_from_base_ring(self):
"""
sage: J(1)
Traceback (most recent call last):
...
- ValueError: not a naturally-represented algebra element
+ ValueError: not an element of this algebra
"""
return None
- def __init__(self,
- field,
- mult_table,
- prefix='e',
- category=None,
- natural_basis=None,
- check=True):
+
+ 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 (
- ....: FiniteDimensionalEuclideanJordanAlgebra,
- ....: 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,z = J.random_elements(3)
+ 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: x*y == y*x
+ 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
+
+ """
+ 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 ``field`` we're given must be real with ``check=True``::
+ The values we've presupplied to the constructors agree with
+ the computation::
- sage: JordanSpinEJA(2,QQbar)
- Traceback (most recent call last):
- ...
- ValueError: field is not real
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: J.is_associative() == J._jordan_product_is_associative()
+ True
- The multiplication table must be square with ``check=True``::
+ """
+ R = self.base_ring()
- sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
- Traceback (most recent call last):
- ...
- ValueError: multiplication table is not square
+ # 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.
"""
- 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._natural_basis = natural_basis
-
- if category is None:
- category = MagmaticAlgebras(field).FiniteDimensional()
- category = category.WithBasis().Unital()
-
- # The multiplication table had better be square
- n = len(mult_table)
- if check:
- if not all( len(l) == n for l in mult_table ):
- raise ValueError("multiplication table is not square")
-
- fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
- fda.__init__(field,
- range(n),
- 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 = [
- list(map(lambda x: self.from_vector(x), ls))
- for ls in mult_table
- ]
-
- if check:
- if not self._is_commutative():
- raise ValueError("algebra is not commutative")
- 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")
+ 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
SETUP::
- sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: JordanSpinEJA,
....: HadamardEJA,
....: RealSymmetricEJA)
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 = HadamardEJA.random_instance()
- sage: x = J.random_element()
- sage: J(x.to_vector().column()) == x
- True
- sage: J = JordanSpinEJA.random_instance()
+ sage: J = random_eja()
sage: x = J.random_element()
- sage: J(x.to_vector().column()) == x
+ sage: J(x.to_matrix()) == x
True
+ We cannot coerce elements between algebras just because their
+ matrix representations are compatible::
+
+ sage: J1 = HadamardEJA(3)
+ sage: J2 = JordanSpinEJA(3)
+ sage: J2(J1.one())
+ Traceback (most recent call last):
+ ...
+ ValueError: not an element of this algebra
+ sage: J1(J2.zero())
+ Traceback (most recent call last):
+ ...
+ ValueError: not an element of this algebra
"""
- msg = "not a naturally-represented algebra element"
- if elt == 0:
- # The superclass implementation of random_element()
- # needs to be able to coerce "0" into the algebra.
- return self.zero()
- elif elt in self.base_ring():
+ 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)
- natural_basis = self.natural_basis()
- basis_space = natural_basis[0].matrix_space()
- if elt not in basis_space:
+ 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
+
+ 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)
- @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.
+ try:
+ coords = W.coordinate_vector(V(elt))
+ except ArithmeticError: # vector is not in free module
+ raise ValueError(msg)
- 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
+ return self.from_vector(coords)
def _repr_(self):
"""
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 _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=False`` and passed
- an invalid multiplication table.
- """
- return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
- for i in range(self.dimension())
- for j in range(self.dimension()) )
-
- 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=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 _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=False`` and passed
- an invalid multiplication table.
- """
-
- # Used to check whether or not something is zero in an inexact
- # ring. This number is sufficient to allow the construction of
- # QuaternionHermitianEJA(2, RDF) with check=True.
- epsilon = 1e-16
-
- 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 self.base_ring().is_exact():
- if diff != 0:
- return False
- else:
- if diff.abs() > epsilon:
- return False
-
- return True
@cached_method
def characteristic_polynomial_of(self):
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::
+
+ 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):
"""
SETUP::
- sage: from mjo.eja.eja_algebra import random_eja
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: HadamardEJA,
+ ....: BilinearFormEJA)
EXAMPLES:
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):
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 range(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, 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.
- 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.)
+ 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.
- Note that this will always return a matrix. The standard basis
- in `R^n` will be returned as `n`-by-`1` column matrices.
+ 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::
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 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):
+ def matrix_space(self):
"""
- Return the matrix space in which this algebra's natural basis
- elements live.
+ 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).
- Generally this will be an `n`-by-`1` column-vector space,
+ "By default" this will be an `n`-by-`1` column-matrix space,
except when the algebra is trivial. There it's `n`-by-`n`
- (where `n` is zero), to ensure that two elements of the
- natural basis space (empty matrices) can be multiplied.
- """
- if self.is_trivial():
- return MatrixSpace(self.base_ring(), 0)
- elif 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()
+ (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.
+ SETUP::
- @staticmethod
- def natural_inner_product(X,Y):
- """
- Compute the inner product of two naturally-represented elements.
+ sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+ ....: JordanSpinEJA,
+ ....: QuaternionHermitianEJA,
+ ....: TrivialEJA)
+
+ EXAMPLES:
+
+ By default, the matrix representation is just a column-matrix
+ equivalent to the vector representation::
+
+ sage: J = JordanSpinEJA(3)
+ sage: J.matrix_space()
+ Full MatrixSpace of 3 by 1 dense matrices over Algebraic
+ Real Field
+
+ The matrix representation in the trivial algebra is
+ zero-by-zero instead of the usual `n`-by-one::
+
+ sage: J = TrivialEJA()
+ sage: J.matrix_space()
+ Full MatrixSpace of 0 by 0 dense matrices over Algebraic
+ Real Field
+
+ The matrix space for complex/quaternion Hermitian matrix EJA
+ is the space in which their real-embeddings live, not the
+ original complex/quaternion matrix space::
+
+ sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+ sage: J.matrix_space()
+ Full MatrixSpace of 4 by 4 dense matrices over Rational Field
+ sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
+ sage: J.matrix_space()
+ Full MatrixSpace of 4 by 4 dense matrices over Rational Field
- 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()
+ if self.is_trivial():
+ return MatrixSpace(self.base_ring(), 0)
+ else:
+ return self.matrix_basis()[0].parent()
@cached_method
sage: from mjo.eja.eja_algebra import (HadamardEJA,
....: random_eja)
- EXAMPLES::
+ EXAMPLES:
+
+ 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()
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
+
+ ::
+
+ 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
# 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):
Vector space of degree 6 and dimension 2...
sage: J1
Euclidean Jordan algebra of dimension 3...
- sage: J0.one().natural_representation()
+ 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]).natural_representation()
+ 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]).natural_representation()
+ 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().natural_representation()
+ sage: J1.one().to_matrix()
[1 0 0]
[0 1 0]
[0 0 0]
# 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 = FiniteDimensionalEuclideanJordanSubalgebra(self, ())
+ trivial = self.subalgebra(())
J0 = trivial # eigenvalue zero
J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
J1 = trivial # eigenvalue one
J5 = eigspace
else:
gens = tuple( self.from_vector(b) for b in eigspace.basis() )
- subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, gens)
+ subalg = self.subalgebra(gens, check_axioms=False)
if eigval == 0:
J0 = subalg
elif eigval == 1:
return tuple( self.random_element(thorough)
for idx in range(count) )
- @classmethod
- def random_instance(cls, field=AA, **kwargs):
- """
- Return a random instance of this type of algebra.
-
- 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)
-
- n = ZZ.random_element(cls._max_test_case_size() + 1)
- return cls(n, field, **kwargs)
@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
+
"""
n = self.dimension()
- var_names = [ "X" + str(z) for z in range(1,n+1) ]
- R = PolynomialRing(self.base_ring(), var_names)
+ R = self.coordinate_polynomial_ring()
vars = R.gens()
F = R.fraction_field()
# The theory says that only the first "r" coefficients are
# nonzero, and they actually live in the original polynomial
- # ring and not the fraction field. We negate them because
- # in the actual characteristic polynomial, they get moved
- # to the other side where x^r lives.
- return -A_rref.solve_right(E*b).change_ring(R)[:r]
+ # ring and not the fraction field. We negate them because in
+ # the actual characteristic polynomial, they get moved to the
+ # other side where x^r lives. We don't bother to trim A_rref
+ # down to a square matrix and solve the resulting system,
+ # because the upper-left r-by-r portion of A_rref is
+ # guaranteed to be the identity matrix, so e.g.
+ #
+ # A_rref.solve_right(Y)
+ #
+ # would just be returning Y.
+ return (-E*b)[:r].change_ring(R)
@cached_method
def rank(self):
Ensure that computing the rank actually works, since the ranks
of all simple algebras are known and will be cached by default::
- sage: J = HadamardEJA(4)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 4
-
- ::
-
- sage: J = JordanSpinEJA(4)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 2
-
- ::
-
- sage: J = RealSymmetricEJA(3)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 3
-
- ::
+ 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
- sage: J = ComplexHermitianEJA(2)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 2
+ """
+ return len(self._charpoly_coefficients())
- ::
- sage: J = QuaternionHermitianEJA(2)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 2
+ def subalgebra(self, basis, **kwargs):
+ r"""
+ Create a subalgebra of this algebra from the given basis.
"""
- return len(self._charpoly_coefficients())
+ from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+ return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
def vector_space(self):
return self.zero().to_vector().parent().ambient_vector_space()
- Element = FiniteDimensionalEuclideanJordanAlgebraElement
-
-class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra):
- """
- Return the Euclidean Jordan Algebra corresponding to the set
- `R^n` under the Hadamard product.
-
- 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.
+class RationalBasisEJA(FiniteDimensionalEJA):
+ r"""
+ New class for algebras whose supplied basis elements have all rational entries.
SETUP::
- sage: from mjo.eja.eja_algebra import HadamardEJA
+ sage: from mjo.eja.eja_algebra import BilinearFormEJA
EXAMPLES:
- This multiplication table can be verified by hand::
-
- sage: J = HadamardEJA(3)
- sage: e0,e1,e2 = J.gens()
- sage: e0*e0
- e0
- sage: e0*e1
- 0
- sage: e0*e2
- 0
- sage: e1*e1
- e1
- sage: e1*e2
- 0
- sage: e2*e2
- e2
-
- TESTS:
-
- We can change the generator prefix::
+ The supplied basis is orthonormalized by default::
- sage: HadamardEJA(3, prefix='r').gens()
- (r0, r1, r2)
+ 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]
+ )
"""
- def __init__(self, n, field=AA, **kwargs):
- V = VectorSpace(field, n)
- mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
- for i in range(n) ]
-
- fdeja = super(HadamardEJA, self)
- fdeja.__init__(field, mult_table, **kwargs)
- self.rank.set_cache(n)
-
- def inner_product(self, x, y):
- """
- Faster to reimplement than to use natural representations.
+ 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.
+ if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
+ raise TypeError("basis not rational")
+
+ super().__init__(basis,
+ jordan_product,
+ inner_product,
+ field=field,
+ check_field=check_field,
+ **kwargs)
+
+ self._rational_algebra = None
+ if field is not QQ:
+ # There's no point in constructing the extra algebra if this
+ # 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,
+ associative=self.is_associative(),
+ orthonormalize=False,
+ check_field=False,
+ check_axioms=False)
+ @cached_method
+ def _charpoly_coefficients(self):
+ r"""
SETUP::
- sage: from mjo.eja.eja_algebra import HadamardEJA
-
- TESTS:
-
- Ensure that this is the usual inner product for the algebras
- over `R^n`::
+ sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
+ ....: JordanSpinEJA)
- sage: set_random_seed()
- sage: J = HadamardEJA.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
+ EXAMPLES:
- """
- return x.to_vector().inner_product(y.to_vector())
+ 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)::
+ 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(RationalBasisEJA):
+ r"""
+ A class for the Euclidean Jordan algebras that we know by name.
-def random_eja(field=AA):
- """
- Return a "random" finite-dimensional Euclidean Jordan Algebra.
+ 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 random_eja
-
- TESTS::
+ sage: from mjo.eja.eja_algebra import ConcreteEJA
- sage: random_eja()
- Euclidean Jordan algebra of dimension...
+ TESTS:
- """
- classname = choice([TrivialEJA,
- HadamardEJA,
- JordanSpinEJA,
- RealSymmetricEJA,
- ComplexHermitianEJA,
- QuaternionHermitianEJA])
- return classname.random_instance(field=field)
+ Our basis is normalized with respect to the algebra's inner
+ product, unless we specify otherwise::
+ sage: set_random_seed()
+ sage: J = ConcreteEJA.random_instance()
+ sage: all( b.norm() == 1 for b in J.gens() )
+ True
+ 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::
+ sage: set_random_seed()
+ sage: J = ConcreteEJA.random_instance()
+ sage: x = J.random_element()
+ sage: x.operator().is_self_adjoint()
+ 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, normalize_basis=True, **kwargs):
+ def _max_random_instance_size():
"""
- 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).
+ 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.
+
+ This method must be implemented in each subclass.
"""
- # Used in this class's fast _charpoly_coefficients() override.
- self._basis_normalizers = None
+ raise NotImplementedError
- # 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)
+ @classmethod
+ def random_instance(cls, *args, **kwargs):
+ """
+ Return a random instance of this type of algebra.
- if len(basis) > 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 zip(basis,self._basis_normalizers))
+ This method should be implemented in each subclass.
+ """
+ from sage.misc.prandom import choice
+ eja_class = choice(cls.__subclasses__())
- Qs = self.multiplication_table_from_matrix_basis(basis)
+ # 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)
- fdeja = super(MatrixEuclideanJordanAlgebra, self)
- fdeja.__init__(field, Qs, natural_basis=basis, **kwargs)
- return
+class MatrixEJA:
+ @staticmethod
+ def jordan_product(X,Y):
+ return (X*Y + Y*X)/2
- @cached_method
- def _charpoly_coefficients(self):
+ @staticmethod
+ def trace_inner_product(X,Y):
r"""
- Override the parent method with something that tries to compute
- over a faster (non-extension) field.
+ A trace inner-product for matrices that aren't embedded in the
+ reals.
"""
- 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_coefficients()
- else:
- basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
- self._basis_normalizers) )
-
- # Do this over the rationals and convert back at the end.
- # Only works because we know the entries of the basis are
- # integers.
- J = MatrixEuclideanJordanAlgebra(QQ,
- basis,
- normalize_basis=False)
- a = J._charpoly_coefficients()
-
- # Unfortunately, changing the basis does change the
- # coefficients of the characteristic polynomial, but since
- # these are really the coefficients of the "characteristic
- # polynomial of" function, everything is still nice and
- # unevaluated. It's therefore "obvious" how scaling the
- # basis affects the coordinate variables X1, X2, et
- # cetera. Scaling the first basis vector up by "n" adds a
- # factor of 1/n into every "X1" term, for example. So here
- # we simply undo the basis_normalizer scaling that we
- # performed earlier.
- #
- # The a[0] access here is safe because trivial algebras
- # won't have any basis normalizers and therefore won't
- # make it to this "else" branch.
- XS = a[0].parent().gens()
- subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i]
- for i in range(len(XS)) }
- return tuple( a_i.subs(subs_dict) for a_i in a )
-
+ # We take the norm (absolute value) because Octonions() isn't
+ # smart enough yet to coerce its one() into the base field.
+ return (X*Y).trace().abs()
+class RealEmbeddedMatrixEJA(MatrixEJA):
@staticmethod
- def 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.
- if len(basis) == 0:
- return []
-
- 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 range(n)] for i in range(n)]
- for i in range(n):
- for j in range(n):
- mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
- mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
+ def dimension_over_reals():
+ r"""
+ The dimension of this matrix's base ring over the reals.
- return mult_table
+ The reals are dimension one over themselves, obviously; that's
+ just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
+ have dimension two. Finally, the quaternions have dimension
+ four over the reals.
+ This is used to determine the size of the matrix returned from
+ :meth:`real_embed`, among other things.
+ """
+ raise NotImplementedError
- @staticmethod
- def real_embed(M):
+ @classmethod
+ def real_embed(cls,M):
"""
Embed the matrix ``M`` into a space of real matrices.
real_embed(M*N) = real_embed(M)*real_embed(N)
"""
- raise NotImplementedError
+ if M.ncols() != M.nrows():
+ raise ValueError("the matrix 'M' must be square")
+ return M
- @staticmethod
- def real_unembed(M):
+ @classmethod
+ def real_unembed(cls,M):
"""
The inverse of :meth:`real_embed`.
"""
- raise NotImplementedError
+ if M.ncols() != M.nrows():
+ raise ValueError("the matrix 'M' must be square")
+ if not ZZ(M.nrows()).mod(cls.dimension_over_reals()).is_zero():
+ raise ValueError("the matrix 'M' must be a real embedding")
+ return M
@classmethod
- def natural_inner_product(cls,X,Y):
- Xu = cls.real_unembed(X)
- Yu = cls.real_unembed(Y)
- tr = (Xu*Yu).trace()
+ def trace_inner_product(cls,X,Y):
+ r"""
+ Compute the trace inner-product of two real-embeddings.
- if tr in RLF:
- # It's real already.
- return tr
+ SETUP::
- # 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.
- return tr.coefficient_tuple()[0]
+ sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+ ....: QuaternionHermitianEJA)
+ EXAMPLES::
-class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
- @staticmethod
- def real_embed(M):
- """
- The identity function, for embedding real matrices into real
- matrices.
- """
- return M
+ sage: set_random_seed()
+ sage: J = ComplexHermitianEJA.random_instance()
+ sage: x,y = J.random_elements(2)
+ sage: Xe = x.to_matrix()
+ sage: Ye = y.to_matrix()
+ sage: X = J.real_unembed(Xe)
+ sage: Y = J.real_unembed(Ye)
+ sage: expected = (X*Y).trace().real()
+ sage: actual = J.trace_inner_product(Xe,Ye)
+ sage: actual == expected
+ True
- @staticmethod
- def real_unembed(M):
- """
- The identity function, for unembedding real matrices from real
- matrices.
- """
- return M
+ ::
+
+ sage: set_random_seed()
+ sage: J = QuaternionHermitianEJA.random_instance()
+ sage: x,y = J.random_elements(2)
+ sage: Xe = x.to_matrix()
+ sage: Ye = y.to_matrix()
+ sage: X = J.real_unembed(Xe)
+ sage: Y = J.real_unembed(Ye)
+ sage: expected = (X*Y).trace().coefficient_tuple()[0]
+ sage: actual = J.trace_inner_product(Xe,Ye)
+ sage: actual == expected
+ True
+ """
+ # This does in fact compute the real part of the trace.
+ # If we compute the trace of e.g. a complex matrix M,
+ # then we do so by adding up its diagonal entries --
+ # call them z_1 through z_n. The real embedding of z_1
+ # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
+ # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
+ return (X*Y).trace()/cls.dimension_over_reals()
-class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
+class RealSymmetricEJA(ConcreteEJA, MatrixEJA):
"""
The rank-n simple EJA consisting of real symmetric n-by-n
matrices, the usual symmetric Jordan product, and the trace inner
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, RDF)
+ sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
Euclidean Jordan algebra of dimension 3 over Real Double Field
- sage: RealSymmetricEJA(2, RR)
+ sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
Euclidean Jordan algebra of dimension 3 over Real Field with
53 bits of precision
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
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
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::
-
- sage: set_random_seed()
- sage: x = RealSymmetricEJA.random_instance().random_element()
- sage: x.operator().matrix().is_symmetric()
- True
-
We can construct the (trivial) algebra of rank zero::
sage: RealSymmetricEJA(0)
sage: set_random_seed()
sage: n = ZZ.random_element(1,5)
- sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
+ sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ)
sage: all( M.is_symmetric() for M in B)
True
else:
Sij = Eij + Eij.transpose()
S.append(Sij)
- return S
+ return tuple(S)
@staticmethod
- def _max_test_case_size():
+ def _max_random_instance_size():
return 4 # Dimension 10
+ @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=AA, **kwargs):
- basis = self._denormalized_basis(n, field)
- super(RealSymmetricEJA, self).__init__(field, basis, **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
+
+ associative = False
+ if n <= 1:
+ associative = True
+
+ super().__init__(self._denormalized_basis(n,field),
+ self.jordan_product,
+ self.trace_inner_product,
+ field=field,
+ associative=associative,
+ **kwargs)
+
+ # TODO: this could be factored out somehow, but is left here
+ # because the MatrixEJA is not presently a subclass of the
+ # FDEJA class that defines rank() and one().
self.rank.set_cache(n)
+ idV = self.matrix_space().one()
+ self.one.set_cache(self(idV))
+
+
+
+class ComplexMatrixEJA(RealEmbeddedMatrixEJA):
+ # A manual dictionary-cache for the complex_extension() method,
+ # since apparently @classmethods can't also be @cached_methods.
+ _complex_extension = {}
+
+ @classmethod
+ def complex_extension(cls,field):
+ r"""
+ The complex field that we embed/unembed, as an extension
+ of the given ``field``.
+ """
+ if field in cls._complex_extension:
+ return cls._complex_extension[field]
+
+ # Sage doesn't know how to adjoin the complex "i" (the root of
+ # x^2 + 1) to a field in a general way. Here, we just enumerate
+ # all of the cases that I have cared to support so far.
+ if field is AA:
+ # Sage doesn't know how to embed AA into QQbar, i.e. how
+ # to adjoin sqrt(-1) to AA.
+ F = QQbar
+ elif not field.is_exact():
+ # RDF or RR
+ F = field.complex_field()
+ else:
+ # Works for QQ and... maybe some other fields.
+ R = PolynomialRing(field, 'z')
+ z = R.gen()
+ F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+ cls._complex_extension[field] = F
+ return F
-class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@staticmethod
- def real_embed(M):
+ def dimension_over_reals():
+ return 2
+
+ @classmethod
+ def real_embed(cls,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 +
SETUP::
- sage: from mjo.eja.eja_algebra import \
- ....: ComplexMatrixEuclideanJordanAlgebra
+ sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
EXAMPLES::
sage: x3 = F(-i)
sage: x4 = F(6)
sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
- sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
+ sage: ComplexMatrixEJA.real_embed(M)
[ 4 -2| 1 2]
[ 2 4|-2 1]
[-----+-----]
Embedding is a homomorphism (isomorphism, in fact)::
sage: set_random_seed()
- sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
- sage: n = ZZ.random_element(n_max)
+ sage: n = ZZ.random_element(3)
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 = ComplexMatrixEJA.real_embed(X)
+ sage: Ye = ComplexMatrixEJA.real_embed(Y)
+ sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
sage: Xe*Ye == XYe
True
"""
+ super().real_embed(M)
n = M.nrows()
- if M.ncols() != n:
- raise ValueError("the matrix 'M' must be square")
# We don't need any adjoined elements...
field = M.base_ring().base_ring()
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]]))
+ a = z.real()
+ b = z.imag()
+ blocks.append(matrix(field, 2, [ [ a, b],
+ [-b, a] ]))
return matrix.block(field, n, blocks)
- @staticmethod
- def real_unembed(M):
+ @classmethod
+ def real_unembed(cls,M):
"""
The inverse of _embed_complex_matrix().
SETUP::
- sage: from mjo.eja.eja_algebra import \
- ....: ComplexMatrixEuclideanJordanAlgebra
+ sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
EXAMPLES::
....: [-2, 1, -4, 3],
....: [ 9, 10, 11, 12],
....: [-10, 9, -12, 11] ])
- sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
+ sage: ComplexMatrixEJA.real_unembed(A)
[ 2*I + 1 4*I + 3]
[ 10*I + 9 12*I + 11]
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
+ sage: Me = ComplexMatrixEJA.real_embed(M)
+ sage: ComplexMatrixEJA.real_unembed(Me) == M
True
"""
+ super().real_unembed(M)
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()
- if field is AA:
- # Sage doesn't know how to embed AA into QQbar, i.e. how
- # to adjoin sqrt(-1) to AA.
- F = QQbar
- else:
- F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+ d = cls.dimension_over_reals()
+ F = cls.complex_extension(M.base_ring())
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 range(n/2):
- for j in range(n/2):
- submat = M[2*k:2*k+2,2*j:2*j+2]
+ for k in range(n/d):
+ for j in range(n/d):
+ submat = M[d*k:d*k+d,d*j:d*j+d]
if submat[0,0] != submat[1,1]:
raise ValueError('bad on-diagonal submatrix')
if submat[0,1] != -submat[1,0]:
z = submat[0,0] + submat[0,1]*i
elements.append(z)
- return matrix(F, n/2, elements)
-
-
- @classmethod
- def natural_inner_product(cls,X,Y):
- """
- 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().real()
- sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
- sage: actual == expected
- True
-
- """
- return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2
+ return matrix(F, n/d, elements)
-class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
+class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
"""
The rank-n simple EJA consisting of complex Hermitian n-by-n
matrices over the real numbers, the usual symmetric Jordan product,
In theory, our "field" can be any subfield of the reals::
- sage: ComplexHermitianEJA(2, RDF)
+ sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
Euclidean Jordan algebra of dimension 4 over Real Double Field
- sage: ComplexHermitianEJA(2, RR)
+ sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
Euclidean Jordan algebra of dimension 4 over Real Field with
53 bits of precision
The dimension of this algebra is `n^2`::
sage: set_random_seed()
- sage: n_max = ComplexHermitianEJA._max_test_case_size()
+ 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
sage: set_random_seed()
sage: J = ComplexHermitianEJA.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
sage: ComplexHermitianEJA(2, prefix='z').gens()
(z0, z1, z2, z3)
- Our natural basis is normalized with respect to the natural inner
- product unless we specify otherwise::
-
- sage: set_random_seed()
- sage: J = ComplexHermitianEJA.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::
-
- sage: set_random_seed()
- sage: x = ComplexHermitianEJA.random_instance().random_element()
- sage: x.operator().matrix().is_symmetric()
- True
-
We can construct the (trivial) algebra of rank zero::
sage: ComplexHermitianEJA(0)
sage: set_random_seed()
sage: n = ZZ.random_element(1,5)
- sage: field = QuadraticField(2, 'sqrt2')
- sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
+ sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ)
sage: all( M.is_symmetric() for M in B)
True
"""
- R = PolynomialRing(field, 'z')
+ R = PolynomialRing(ZZ, 'z')
z = R.gen()
- F = field.extension(z**2 + 1, 'I')
- I = F.gen()
+ F = ZZ.extension(z**2 + 1, 'I')
+ I = F.gen(1)
# This is like the symmetric case, but we need to be careful:
#
# * The diagonal will (as a result) be real.
#
S = []
+ Eij = matrix.zero(F,n)
for i in range(n):
for j in range(i+1):
- Eij = matrix(F, n, lambda k,l: k==i and l==j)
+ # "build" E_ij
+ Eij[i,j] = 1
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())
+ Eij[j,i] = 1 # Eij = Eij + Eij.transpose()
+ Sij_real = cls.real_embed(Eij)
S.append(Sij_real)
- Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
+ # Eij = I*Eij - I*Eij.transpose()
+ Eij[i,j] = I
+ Eij[j,i] = -I
+ Sij_imag = cls.real_embed(Eij)
S.append(Sij_imag)
+ Eij[j,i] = 0
+ # "erase" E_ij
+ Eij[i,j] = 0
- # Since we embedded these, we can drop back to the "field" that we
- # started with instead of the complex extension "F".
- return ( s.change_ring(field) for s in S )
+ # Since we embedded the entries, we can drop back to the
+ # desired real "field" instead of the extension "F".
+ return tuple( s.change_ring(field) for s in S )
def __init__(self, n, field=AA, **kwargs):
- basis = self._denormalized_basis(n,field)
- super(ComplexHermitianEJA,self).__init__(field, basis, **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
+
+ associative = False
+ if n <= 1:
+ associative = True
+
+ super().__init__(self._denormalized_basis(n,field),
+ self.jordan_product,
+ self.trace_inner_product,
+ field=field,
+ associative=associative,
+ **kwargs)
+ # TODO: this could be factored out somehow, but is left here
+ # because the MatrixEJA is not presently a subclass of the
+ # FDEJA class that defines rank() and one().
self.rank.set_cache(n)
+ idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+ self.one.set_cache(self(idV))
+
+ @staticmethod
+ def _max_random_instance_size():
+ return 3 # Dimension 9
+
+ @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)
+
+class QuaternionMatrixEJA(RealEmbeddedMatrixEJA):
+ # A manual dictionary-cache for the quaternion_extension() method,
+ # since apparently @classmethods can't also be @cached_methods.
+ _quaternion_extension = {}
+
+ @classmethod
+ def quaternion_extension(cls,field):
+ r"""
+ The quaternion field that we embed/unembed, as an extension
+ of the given ``field``.
+ """
+ if field in cls._quaternion_extension:
+ return cls._quaternion_extension[field]
+
+ Q = QuaternionAlgebra(field,-1,-1)
+
+ cls._quaternion_extension[field] = Q
+ return Q
-class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@staticmethod
- def real_embed(M):
+ def dimension_over_reals():
+ return 4
+
+ @classmethod
+ def real_embed(cls,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
SETUP::
- sage: from mjo.eja.eja_algebra import \
- ....: QuaternionMatrixEuclideanJordanAlgebra
+ sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
EXAMPLES::
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)
+ sage: QuaternionMatrixEJA.real_embed(M)
[ 1 2 3 4]
[-2 1 -4 3]
[-3 4 1 -2]
Embedding is a homomorphism (isomorphism, in fact)::
sage: set_random_seed()
- sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
- sage: n = ZZ.random_element(n_max)
+ sage: n = ZZ.random_element(2)
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 = QuaternionMatrixEJA.real_embed(X)
+ sage: Ye = QuaternionMatrixEJA.real_embed(Y)
+ sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
sage: Xe*Ye == XYe
True
"""
+ super().real_embed(M)
quaternions = M.base_ring()
n = M.nrows()
- if M.ncols() != n:
- raise ValueError("the matrix 'M' must be square")
F = QuadraticField(-1, 'I')
i = F.gen()
d = t[3]
cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
[-c + d*i, a - b*i]])
- realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM)
+ realM = ComplexMatrixEJA.real_embed(cplxM)
blocks.append(realM)
# We should have real entries by now, so use the realest field
- @staticmethod
- def real_unembed(M):
+ @classmethod
+ def real_unembed(cls,M):
"""
The inverse of _embed_quaternion_matrix().
SETUP::
- sage: from mjo.eja.eja_algebra import \
- ....: QuaternionMatrixEuclideanJordanAlgebra
+ sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
EXAMPLES::
....: [-2, 1, -4, 3],
....: [-3, 4, 1, -2],
....: [-4, -3, 2, 1]])
- sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
+ sage: QuaternionMatrixEJA.real_unembed(M)
[1 + 2*i + 3*j + 4*k]
TESTS:
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
+ sage: Me = QuaternionMatrixEJA.real_embed(M)
+ sage: QuaternionMatrixEJA.real_unembed(Me) == M
True
"""
+ super().real_unembed(M)
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")
+ d = cls.dimension_over_reals()
# 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 range(n/4):
- for m in range(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].real()
- z += submat[0,0].imag()*i
- z += submat[0,1].real()*j
- z += submat[0,1].imag()*k
- elements.append(z)
-
- return matrix(Q, n/4, elements)
-
-
- @classmethod
- def natural_inner_product(cls,X,Y):
- """
- 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
+ Q = cls.quaternion_extension(M.base_ring())
+ i,j,k = Q.gens()
- """
- return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4
+ # 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 range(n/d):
+ for m in range(n/d):
+ submat = ComplexMatrixEJA.real_unembed(
+ M[d*l:d*l+d,d*m:d*m+d] )
+ 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].real()
+ z += submat[0,0].imag()*i
+ z += submat[0,1].real()*j
+ z += submat[0,1].imag()*k
+ elements.append(z)
+ return matrix(Q, n/d, elements)
-class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
- """
+
+class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
+ 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
In theory, our "field" can be any subfield of the reals::
- sage: QuaternionHermitianEJA(2, RDF)
+ sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
Euclidean Jordan algebra of dimension 6 over Real Double Field
- sage: QuaternionHermitianEJA(2, RR)
+ sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
Euclidean Jordan algebra of dimension 6 over Real Field with
53 bits of precision
The dimension of this algebra is `2*n^2 - n`::
sage: set_random_seed()
- sage: n_max = QuaternionHermitianEJA._max_test_case_size()
+ sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
sage: n = ZZ.random_element(1, n_max)
sage: J = QuaternionHermitianEJA(n)
sage: J.dimension() == 2*(n^2) - n
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: 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: 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::
-
- sage: set_random_seed()
- sage: J = QuaternionHermitianEJA.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::
-
- sage: set_random_seed()
- sage: x = QuaternionHermitianEJA.random_instance().random_element()
- sage: x.operator().matrix().is_symmetric()
- True
-
We can construct the (trivial) algebra of rank zero::
sage: QuaternionHermitianEJA(0)
sage: set_random_seed()
sage: n = ZZ.random_element(1,5)
- sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
+ sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ)
sage: all( M.is_symmetric() for M in B )
True
# * The diagonal will (as a result) be real.
#
S = []
+ Eij = matrix.zero(Q,n)
for i in range(n):
for j in range(i+1):
- Eij = matrix(Q, n, lambda k,l: k==i and l==j)
+ # "build" E_ij
+ Eij[i,j] = 1
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())
+ # Eij = Eij + Eij.transpose()
+ Eij[j,i] = 1
+ Sij_real = cls.real_embed(Eij)
S.append(Sij_real)
- Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
+ # Eij = I*(Eij - Eij.transpose())
+ Eij[i,j] = I
+ Eij[j,i] = -I
+ Sij_I = cls.real_embed(Eij)
S.append(Sij_I)
- Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
+ # Eij = J*(Eij - Eij.transpose())
+ Eij[i,j] = J
+ Eij[j,i] = -J
+ Sij_J = cls.real_embed(Eij)
S.append(Sij_J)
- Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
+ # Eij = K*(Eij - Eij.transpose())
+ Eij[i,j] = K
+ Eij[j,i] = -K
+ Sij_K = cls.real_embed(Eij)
S.append(Sij_K)
+ Eij[j,i] = 0
+ # "erase" E_ij
+ Eij[i,j] = 0
+
+ # Since we embedded the entries, we can drop back to the
+ # desired real "field" instead of the quaternion algebra "Q".
+ return tuple( s.change_ring(field) for s in S )
+
+
+ def __init__(self, n, 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
+
+ associative = False
+ if n <= 1:
+ associative = True
+
+ super().__init__(self._denormalized_basis(n,field),
+ self.jordan_product,
+ self.trace_inner_product,
+ field=field,
+ associative=associative,
+ **kwargs)
+
+ # TODO: this could be factored out somehow, but is left here
+ # because the MatrixEJA is not presently a subclass of the
+ # FDEJA class that defines rank() and one().
+ self.rank.set_cache(n)
+ idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+ self.one.set_cache(self(idV))
+
+
+ @staticmethod
+ def _max_random_instance_size():
+ r"""
+ The maximum rank of a random QuaternionHermitianEJA.
+ """
+ return 2 # Dimension 6
+
+ @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)
+
+
+class HadamardEJA(ConcreteEJA):
+ """
+ Return the Euclidean Jordan Algebra corresponding to the set
+ `R^n` under the Hadamard product.
+
+ 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.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import HadamardEJA
+
+ EXAMPLES:
+
+ This multiplication table can be verified by hand::
+
+ 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
+
+ TESTS:
- # 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 )
+ We can change the generator prefix::
+ sage: HadamardEJA(3, prefix='r').gens()
+ (r0, r1, r2)
+ """
def __init__(self, n, field=AA, **kwargs):
- basis = self._denormalized_basis(n,field)
- super(QuaternionHermitianEJA,self).__init__(field, basis, **kwargs)
+ if n == 0:
+ jordan_product = lambda x,y: x
+ inner_product = lambda x,y: x
+ else:
+ def jordan_product(x,y):
+ P = x.parent()
+ return P( 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( b.column()
+ for b in FreeModule(field, n).basis() )
+ super().__init__(column_basis,
+ jordan_product,
+ inner_product,
+ field=field,
+ associative=True,
+ **kwargs)
self.rank.set_cache(n)
+ if n == 0:
+ self.one.set_cache( self.zero() )
+ else:
+ self.one.set_cache( sum(self.gens()) )
+
+ @staticmethod
+ def _max_random_instance_size():
+ r"""
+ The maximum dimension of a random HadamardEJA.
+ """
+ return 5
+
+ @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)
-class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
+
+class BilinearFormEJA(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 =
- (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
- symmetric positive-definite "bilinear form" matrix. It has
- dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
- when ``B`` is the identity matrix of order ``n-1``.
+ (<Bx,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::
When no bilinear form is specified, the identity matrix is used,
and the resulting algebra is the Jordan spin algebra::
- sage: J0 = BilinearFormEJA(3)
+ sage: B = matrix.identity(AA,3)
+ sage: J0 = BilinearFormEJA(B)
sage: J1 = JordanSpinEJA(3)
sage: J0.multiplication_table() == J0.multiplication_table()
True
+ An error is raised if the matrix `B` does not correspond to a
+ positive-definite bilinear form::
+
+ 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
+
TESTS:
We can create a zero-dimensional algebra::
- sage: J = BilinearFormEJA(0)
+ sage: B = matrix.identity(AA,0)
+ sage: J = BilinearFormEJA(B)
sage: J.basis()
Finite family {}
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::
+ 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: n = ZZ.random_element(5)
sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
- sage: B = M.transpose()*M
- sage: J = BilinearFormEJA(n, B=B)
+ 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.from_vector(V([0] + (M.inverse()*ei).list()))
+ 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: actual == expected
True
- """
- def __init__(self, n, field=AA, B=None, **kwargs):
- if B is None:
- self._B = matrix.identity(field, max(0,n-1))
- else:
- self._B = B
- V = VectorSpace(field, n)
- mult_table = [[V.zero() for j in range(n)] for i in range(n)]
- for i in range(n):
- for j in range(n):
- x = V.gen(i)
- y = V.gen(j)
- x0 = x[0]
- xbar = x[1:]
- y0 = y[0]
- ybar = y[1:]
- z0 = x0*y0 + (self._B*xbar).inner_product(ybar)
- zbar = y0*xbar + x0*ybar
- z = V([z0] + zbar.list())
- mult_table[i][j] = z
+ """
+ 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
+
+ def inner_product(x,y):
+ return (y.T*B*x)[0,0]
+
+ def jordan_product(x,y):
+ P = x.parent()
+ 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 P([z0] + zbar.list())
+
+ n = B.nrows()
+ column_basis = tuple( b.column()
+ for b in FreeModule(field, n).basis() )
+
+ # TODO: I haven't actually checked this, but it seems legit.
+ associative = False
+ if n <= 2:
+ associative = True
+
+ super().__init__(column_basis,
+ jordan_product,
+ inner_product,
+ field=field,
+ associative=associative,
+ **kwargs)
# The rank of this algebra is two, unless we're in a
# one-dimensional ambient space (because the rank is bounded
# by the ambient dimension).
- fdeja = super(BilinearFormEJA, self)
- fdeja.__init__(field, mult_table, **kwargs)
self.rank.set_cache(min(n,2))
- def inner_product(self, x, y):
- r"""
- Half of the trace inner product.
-
- This is defined so that the special case of the Jordan spin
- algebra gets the usual inner product.
-
- SETUP::
+ if n == 0:
+ self.one.set_cache( self.zero() )
+ else:
+ self.one.set_cache( self.monomial(0) )
- sage: from mjo.eja.eja_algebra import BilinearFormEJA
+ @staticmethod
+ def _max_random_instance_size():
+ r"""
+ The maximum dimension of a random BilinearFormEJA.
+ """
+ return 5
- TESTS:
+ @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)
- Ensure that this is one-half of the trace inner-product when
- the algebra isn't just the reals (when ``n`` isn't one). This
- is in Faraut and Koranyi, and also my "On the symmetry..."
- paper::
+ 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
- sage: set_random_seed()
- sage: n = ZZ.random_element(2,5)
- sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
- sage: B = M.transpose()*M
- sage: J = BilinearFormEJA(n, B=B)
- sage: x = J.random_element()
- sage: y = J.random_element()
- sage: x.inner_product(y) == (x*y).trace()/2
- True
+ from sage.matrix.special import block_matrix
+ B = block_matrix(2,2, [ [B11, ZZ(0) ],
+ [ZZ(0), B22 ] ])
- """
- xvec = x.to_vector()
- xbar = xvec[1:]
- yvec = y.to_vector()
- ybar = yvec[1:]
- return x[0]*y[0] + (self._B*xbar).inner_product(ybar)
+ return cls(B, **kwargs)
class JordanSpinEJA(BilinearFormEJA):
This multiplication table can be verified by hand::
sage: J = JordanSpinEJA(4)
- sage: e0,e1,e2,e3 = J.gens()
- sage: e0*e0
- e0
- sage: e0*e1
- e1
- sage: e0*e2
- e2
- sage: e0*e3
- e3
- sage: e1*e2
+ sage: b0,b1,b2,b3 = J.gens()
+ sage: b0*b0
+ b0
+ sage: b0*b1
+ b1
+ sage: b0*b2
+ b2
+ sage: b0*b3
+ b3
+ sage: b1*b2
0
- sage: e1*e3
+ sage: b1*b3
0
- sage: e2*e3
+ sage: b2*b3
0
We can change the generator prefix::
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, field=AA, **kwargs):
- # This is a special case of the BilinearFormEJA with the identity
- # matrix as its bilinear form.
- return super(JordanSpinEJA, self).__init__(n, field, **kwargs)
+ def __init__(self, n, *args, **kwargs):
+ # This is a special case of the BilinearFormEJA with the
+ # identity matrix as its bilinear form.
+ B = matrix.identity(ZZ, n)
+
+ # 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 5
+
+ @classmethod
+ def random_instance(cls, **kwargs):
+ """
+ Return a random instance of this type of algebra.
+
+ Needed here to override the implementation for ``BilinearFormEJA``.
+ """
+ n = ZZ.random_element(cls._max_random_instance_size() + 1)
+ return cls(n, **kwargs)
-class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
+class TrivialEJA(ConcreteEJA):
"""
The trivial Euclidean Jordan algebra consisting of only a zero element.
0
"""
- def __init__(self, field=AA, **kwargs):
- mult_table = []
- fdeja = super(TrivialEJA, self)
+ def __init__(self, **kwargs):
+ jordan_product = lambda x,y: x
+ inner_product = lambda x,y: 0
+ basis = ()
+
+ # 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,
+ **kwargs)
+
# The rank is zero using my definition, namely the dimension of the
# largest subalgebra generated by any element.
- fdeja.__init__(field, mult_table, **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 DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
+class CartesianProductEJA(FiniteDimensionalEJA):
r"""
- The external (orthogonal) direct sum of two other Euclidean Jordan
- algebras. Essentially the Cartesian product of its two factors.
- Every Euclidean Jordan algebra decomposes into an orthogonal
- direct sum of simple Euclidean Jordan algebras, so no generality
- is lost by providing only this construction.
+ The external (orthogonal) direct sum of two or more Euclidean
+ Jordan algebras. Every Euclidean Jordan algebra decomposes into an
+ orthogonal direct sum of simple Euclidean Jordan algebras which is
+ then isometric to a Cartesian product, so no generality is lost by
+ providing only this construction.
SETUP::
- sage: from mjo.eja.eja_algebra import (HadamardEJA,
- ....: RealSymmetricEJA,
- ....: DirectSumEJA)
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: CartesianProductEJA,
+ ....: HadamardEJA,
+ ....: JordanSpinEJA,
+ ....: RealSymmetricEJA)
- EXAMPLES::
+ EXAMPLES:
+
+ The Jordan product is inherited from our factors and implemented by
+ our CombinatorialFreeModule Cartesian product superclass::
+
+ sage: set_random_seed()
+ sage: J1 = HadamardEJA(2)
+ sage: J2 = RealSymmetricEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: x,y = J.random_elements(2)
+ sage: x*y in J
+ True
+
+ The ability to retrieve the original factors is implemented by our
+ CombinatorialFreeModule Cartesian product superclass::
+
+ sage: J1 = HadamardEJA(2, field=QQ)
+ sage: J2 = JordanSpinEJA(3, field=QQ)
+ sage: J = cartesian_product([J1,J2])
+ sage: J.cartesian_factors()
+ (Euclidean Jordan algebra of dimension 2 over Rational Field,
+ Euclidean Jordan algebra of dimension 3 over Rational Field)
+
+ You can provide more than two factors::
+
+ sage: J1 = HadamardEJA(2)
+ sage: J2 = 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 )
+
+ MS = self.matrix_space()
+ basis = []
+ zero = MS.zero()
+ for i in range(m):
+ for b in factors[i].matrix_basis():
+ z = list(zero)
+ z[i] = b
+ basis.append(z)
+
+ basis = tuple( MS(b) for b in basis )
+
+ # Define jordan/inner products that operate on that matrix_basis.
+ def jordan_product(x,y):
+ return MS(tuple(
+ (factors[i](x[i])*factors[i](y[i])).to_matrix()
+ for i in range(m)
+ ))
+
+ def inner_product(x, y):
+ return sum(
+ factors[i](x[i]).inner_product(factors[i](y[i]))
+ for i in range(m)
+ )
+
+ # There's no need to check the field since it already came
+ # from an EJA. Likewise the axioms are guaranteed to be
+ # satisfied, unless the guy writing this class sucks.
+ #
+ # If you want the basis to be orthonormalized, orthonormalize
+ # the factors.
+ FiniteDimensionalEJA.__init__(self,
+ basis,
+ jordan_product,
+ inner_product,
+ field=field,
+ orthonormalize=False,
+ associative=associative,
+ cartesian_product=True,
+ check_field=False,
+ check_axioms=False)
+
+ ones = tuple(J.one().to_matrix() for J in factors)
+ self.one.set_cache(self(ones))
+ self.rank.set_cache(sum(J.rank() for J in factors))
+
+ def cartesian_factors(self):
+ # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+ return self._sets
+
+ def cartesian_factor(self, i):
+ r"""
+ Return the ``i``th factor of this algebra.
+ """
+ return self._sets[i]
+
+ def _repr_(self):
+ # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+ from sage.categories.cartesian_product import cartesian_product
+ return cartesian_product.symbol.join("%s" % factor
+ for factor in self._sets)
+
+ def matrix_space(self):
+ r"""
+ Return the space that our matrix basis lives in as a Cartesian
+ product.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (HadamardEJA,
+ ....: 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)
+
+ """
+ from sage.categories.cartesian_product import cartesian_product
+ return cartesian_product( [J.matrix_space()
+ for J in self.cartesian_factors()] )
+
+ @cached_method
+ def cartesian_projection(self, i):
+ r"""
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: JordanSpinEJA,
+ ....: HadamardEJA,
+ ....: RealSymmetricEJA,
+ ....: ComplexHermitianEJA)
+
+ EXAMPLES:
+
+ The projection morphisms are Euclidean Jordan algebra
+ operators::
+
+ sage: J1 = HadamardEJA(2)
+ sage: J2 = RealSymmetricEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: J.cartesian_projection(0)
+ Linear operator between finite-dimensional Euclidean Jordan
+ algebras represented by the matrix:
+ [1 0 0 0 0]
+ [0 1 0 0 0]
+ Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
+ Real Field (+) Euclidean Jordan algebra of dimension 3 over
+ Algebraic Real Field
+ Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
+ Real Field
+ sage: J.cartesian_projection(1)
+ Linear operator between finite-dimensional Euclidean Jordan
+ algebras represented by the matrix:
+ [0 0 1 0 0]
+ [0 0 0 1 0]
+ [0 0 0 0 1]
+ Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
+ Real Field (+) Euclidean Jordan algebra of dimension 3 over
+ Algebraic Real Field
+ Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
+ Real Field
+
+ The projections work the way you'd expect on the vector
+ representation of an element::
+
+ sage: J1 = JordanSpinEJA(2)
+ sage: J2 = ComplexHermitianEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: pi_left = J.cartesian_projection(0)
+ sage: pi_right = J.cartesian_projection(1)
+ sage: pi_left(J.one()).to_vector()
+ (1, 0)
+ sage: pi_right(J.one()).to_vector()
+ (1, 0, 0, 1)
+ sage: J.one().to_vector()
+ (1, 0, 1, 0, 0, 1)
+
+ TESTS:
+
+ The answer never changes::
+
+ sage: set_random_seed()
+ sage: J1 = random_eja()
+ sage: J2 = random_eja()
+ sage: J = cartesian_product([J1,J2])
+ sage: P0 = J.cartesian_projection(0)
+ sage: P1 = J.cartesian_projection(0)
+ sage: P0 == P1
+ True
+
+ """
+ offset = sum( self.cartesian_factor(k).dimension()
+ for k in range(i) )
+ Ji = self.cartesian_factor(i)
+ Pi = self._module_morphism(lambda j: Ji.monomial(j - offset),
+ codomain=Ji)
+
+ return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
+
+ @cached_method
+ def cartesian_embedding(self, i):
+ r"""
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: JordanSpinEJA,
+ ....: HadamardEJA,
+ ....: RealSymmetricEJA)
+
+ EXAMPLES:
+
+ The embedding morphisms are Euclidean Jordan algebra
+ operators::
+
+ sage: J1 = HadamardEJA(2)
+ sage: J2 = RealSymmetricEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: J.cartesian_embedding(0)
+ Linear operator between finite-dimensional Euclidean Jordan
+ algebras represented by the matrix:
+ [1 0]
+ [0 1]
+ [0 0]
+ [0 0]
+ [0 0]
+ Domain: Euclidean Jordan algebra of dimension 2 over
+ Algebraic Real Field
+ Codomain: Euclidean Jordan algebra of dimension 2 over
+ Algebraic Real Field (+) Euclidean Jordan algebra of
+ dimension 3 over Algebraic Real Field
+ sage: J.cartesian_embedding(1)
+ Linear operator between finite-dimensional Euclidean Jordan
+ algebras represented by the matrix:
+ [0 0 0]
+ [0 0 0]
+ [1 0 0]
+ [0 1 0]
+ [0 0 1]
+ Domain: Euclidean Jordan algebra of dimension 3 over
+ Algebraic Real Field
+ Codomain: Euclidean Jordan algebra of dimension 2 over
+ Algebraic Real Field (+) Euclidean Jordan algebra of
+ dimension 3 over Algebraic Real Field
+
+ The embeddings work the way you'd expect on the vector
+ representation of an element::
+
+ sage: J1 = JordanSpinEJA(3)
+ sage: J2 = RealSymmetricEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: iota_left = J.cartesian_embedding(0)
+ sage: iota_right = J.cartesian_embedding(1)
+ sage: iota_left(J1.zero()) == J.zero()
+ True
+ sage: iota_right(J2.zero()) == J.zero()
+ True
+ sage: J1.one().to_vector()
+ (1, 0, 0)
+ sage: iota_left(J1.one()).to_vector()
+ (1, 0, 0, 0, 0, 0)
+ sage: J2.one().to_vector()
+ (1, 0, 1)
+ sage: iota_right(J2.one()).to_vector()
+ (0, 0, 0, 1, 0, 1)
+ sage: J.one().to_vector()
+ (1, 0, 0, 1, 0, 1)
+
+ TESTS:
+
+ The answer never changes::
+
+ sage: set_random_seed()
+ sage: J1 = random_eja()
+ sage: J2 = random_eja()
+ sage: J = cartesian_product([J1,J2])
+ sage: E0 = J.cartesian_embedding(0)
+ sage: E1 = J.cartesian_embedding(0)
+ sage: E0 == E1
+ True
+
+ Composing a projection with the corresponding inclusion should
+ produce the identity map, and mismatching them should produce
+ the zero map::
+
+ sage: set_random_seed()
+ sage: J1 = random_eja()
+ sage: J2 = random_eja()
+ sage: J = cartesian_product([J1,J2])
+ sage: iota_left = J.cartesian_embedding(0)
+ sage: iota_right = J.cartesian_embedding(1)
+ sage: pi_left = J.cartesian_projection(0)
+ sage: pi_right = J.cartesian_projection(1)
+ sage: pi_left*iota_left == J1.one().operator()
+ True
+ sage: pi_right*iota_right == J2.one().operator()
+ True
+ sage: (pi_left*iota_right).is_zero()
+ True
+ sage: (pi_right*iota_left).is_zero()
+ True
+
+ """
+ offset = sum( self.cartesian_factor(k).dimension()
+ for k in range(i) )
+ Ji = self.cartesian_factor(i)
+ Ei = Ji._module_morphism(lambda j: self.monomial(j + offset),
+ codomain=self)
+ return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
+
+
+
+FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
+
+class RationalBasisCartesianProductEJA(CartesianProductEJA,
+ RationalBasisEJA):
+ r"""
+ A separate class for products of algebras for which we know a
+ rational basis.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+ ....: RealSymmetricEJA)
+
+ EXAMPLES:
+
+ This gives us fast characteristic polynomial computations in
+ product algebras, too::
+
+
+ sage: J1 = JordanSpinEJA(2)
sage: J2 = RealSymmetricEJA(3)
- sage: J = DirectSumEJA(J1,J2)
- sage: J.dimension()
- 8
+ sage: J = cartesian_product([J1,J2])
+ sage: J.characteristic_polynomial_of().degree()
+ 5
sage: J.rank()
5
"""
- def __init__(self, J1, J2, field=AA, **kwargs):
- n1 = J1.dimension()
- n2 = J2.dimension()
- n = n1+n2
- V = VectorSpace(field, n)
- mult_table = [ [ V.zero() for j in range(n) ]
- for i in range(n) ]
- for i in range(n1):
- for j in range(n1):
- p = (J1.monomial(i)*J1.monomial(j)).to_vector()
- mult_table[i][j] = V(p.list() + [field.zero()]*n2)
-
- for i in range(n2):
- for j in range(n2):
- p = (J2.monomial(i)*J2.monomial(j)).to_vector()
- mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
-
- fdeja = super(DirectSumEJA, self)
- fdeja.__init__(field, mult_table, **kwargs)
- self.rank.set_cache(J1.rank() + J2.rank())
+ def __init__(self, algebras, **kwargs):
+ CartesianProductEJA.__init__(self, algebras, **kwargs)
+
+ self._rational_algebra = None
+ if self.vector_space().base_field() is not QQ:
+ self._rational_algebra = cartesian_product([
+ r._rational_algebra for r in algebras
+ ])
+
+
+RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
+
+def random_eja(*args, **kwargs):
+ J1 = ConcreteEJA.random_instance(*args, **kwargs)
+
+ # This might make Cartesian products appear roughly as often as
+ # any other ConcreteEJA.
+ if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0:
+ # Use random_eja() again so we can get more than two factors.
+ J2 = random_eja(*args, **kwargs)
+ J = cartesian_product([J1,J2])
+ return J
+ else:
+ return J1
projects / sage.d.git / blobdiff