"""
-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`
+ * :class:`OctonionHermitianEJA`
+
+In addition to these, we provide two other example constructions,
+
+ * :class:`JordanSpinEJA`
+ * :class:`HadamardEJA`
+ * :class:`AlbertEJA`
+ * :class:`TrivialEJA`
+
+The Jordan spin algebra is a bilinear form algebra where the bilinear
+form is the identity. The Hadamard EJA is simply a Cartesian product
+of one-dimensional spin algebras. The Albert EJA is simply a special
+case of the :class:`OctonionHermitianEJA` where the matrices are
+three-by-three and the resulting space has dimension 27. And
+last/least, the trivial EJA is exactly what you think it is; it could
+also be obtained by constructing a dimension-zero instance of any of
+the other algebras. Cartesian products of these are also supported
+using the usual ``cartesian_product()`` function; as a result, we
+support (up to isomorphism) all Euclidean Jordan algebras.
+
+At a minimum, the following are required to construct a Euclidean
+Jordan algebra:
+
+ * A basis of matrices, column vectors, or MatrixAlgebra elements
+ * A Jordan product defined on the basis
+ * Its inner product defined on the basis
+
+The real numbers form a Euclidean Jordan algebra when both the Jordan
+and inner products are the usual multiplication. We use this as our
+example, and demonstrate a few ways to construct an EJA.
+
+First, we can use one-by-one SageMath matrices with algebraic real
+entries to represent real numbers. We define the Jordan and inner
+products to be essentially real-number multiplication, with the only
+difference being that the Jordan product again returns a one-by-one
+matrix, whereas the inner product must return a scalar. Our basis for
+the one-by-one matrices is of course the set consisting of a single
+matrix with its sole entry non-zero::
+
+ sage: from mjo.eja.eja_algebra import FiniteDimensionalEJA
+ sage: jp = lambda X,Y: X*Y
+ sage: ip = lambda X,Y: X[0,0]*Y[0,0]
+ sage: b1 = matrix(AA, [[1]])
+ sage: J1 = FiniteDimensionalEJA((b1,), jp, ip)
+ sage: J1
+ Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+
+In fact, any positive scalar multiple of that inner-product would work::
+
+ sage: ip2 = lambda X,Y: 16*ip(X,Y)
+ sage: J2 = FiniteDimensionalEJA((b1,), jp, ip2)
+ sage: J2
+ Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+
+But beware that your basis will be orthonormalized _with respect to the
+given inner-product_ unless you pass ``orthonormalize=False`` to the
+constructor. For example::
+
+ sage: J3 = FiniteDimensionalEJA((b1,), jp, ip2, orthonormalize=False)
+ sage: J3
+ Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+
+To see the difference, you can take the first and only basis element
+of the resulting algebra, and ask for it to be converted back into
+matrix form::
+
+ sage: J1.basis()[0].to_matrix()
+ [1]
+ sage: J2.basis()[0].to_matrix()
+ [1/4]
+ sage: J3.basis()[0].to_matrix()
+ [1]
+
+Since square roots are used in that process, the default scalar field
+that we use is the field of algebraic real numbers, ``AA``. You can
+also Use rational numbers, but only if you either pass
+``orthonormalize=False`` or know that orthonormalizing your basis
+won't stray beyond the rational numbers. The example above would
+have worked only because ``sqrt(16) == 4`` is rational.
+
+Another option for your basis is to use elemebts of a
+:class:`MatrixAlgebra`::
+
+ sage: from mjo.matrix_algebra import MatrixAlgebra
+ sage: A = MatrixAlgebra(1,AA,AA)
+ sage: J4 = FiniteDimensionalEJA(A.gens(), jp, ip)
+ sage: J4
+ Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+ sage: J4.basis()[0].to_matrix()
+ +---+
+ | 1 |
+ +---+
+
+An easier way to view the entire EJA basis in its original (but
+perhaps orthonormalized) matrix form is to use the ``matrix_basis``
+method::
+
+ sage: J4.matrix_basis()
+ (+---+
+ | 1 |
+ +---+,)
+
+In particular, a :class:`MatrixAlgebra` is needed to work around the
+fact that matrices in SageMath must have entries in the same
+(commutative and associative) ring as its scalars. There are many
+Euclidean Jordan algebras whose elements are matrices that violate
+those assumptions. The complex, quaternion, and octonion Hermitian
+matrices all have entries in a ring (the complex numbers, quaternions,
+or octonions...) that differs from the algebra's scalar ring (the real
+numbers). Quaternions are also non-commutative; the octonions are
+neither commutative nor associative.
SETUP::
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.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
PolynomialRing,
QuadraticField)
-from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
-from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
-from mjo.eja.eja_utils import _mat2vec
+from mjo.eja.eja_element import (CartesianProductEJAElement,
+ FiniteDimensionalEJAElement)
+from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
+from mjo.eja.eja_utils import _all2list
+
+def EuclideanJordanAlgebras(field):
+ r"""
+ The category of Euclidean Jordan algebras over ``field``, which
+ must be a subfield of the real numbers. For now this is just a
+ convenient wrapper around all of the other category axioms that
+ apply to all EJAs.
+ """
+ category = MagmaticAlgebras(field).FiniteDimensional()
+ category = category.WithBasis().Unital().Commutative()
+ return category
-class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
+class FiniteDimensionalEJA(CombinatorialFreeModule):
r"""
- The lowest-level class for representing a Euclidean Jordan algebra.
+ A finite-dimensional Euclidean Jordan algebra.
+
+ INPUT:
+
+ - ``basis`` -- a tuple; a tuple of basis elements in "matrix
+ form," which must be the same form as the arguments to
+ ``jordan_product`` and ``inner_product``. In reality, "matrix
+ form" can be either vectors, matrices, or a Cartesian product
+ (ordered tuple) of vectors or matrices. All of these would
+ ideally be vector spaces in sage with no special-casing
+ needed; but in reality we turn vectors into column-matrices
+ and Cartesian products `(a,b)` into column matrices
+ `(a,b)^{T}` after converting `a` and `b` themselves.
+
+ - ``jordan_product`` -- a function; afunction of two ``basis``
+ elements (in matrix form) that returns their jordan product,
+ also in matrix form; this will be applied to ``basis`` to
+ compute a multiplication table for the algebra.
+
+ - ``inner_product`` -- a function; a function of two ``basis``
+ elements (in matrix form) that returns their inner
+ product. This will be applied to ``basis`` to compute an
+ inner-product table (basically a matrix) for this algebra.
+
+ - ``matrix_space`` -- the space that your matrix basis lives in,
+ or ``None`` (the default). So long as your basis does not have
+ length zero you can omit this. But in trivial algebras, it is
+ required.
+
+ - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
+ field for the algebra.
+
+ - ``orthonormalize`` -- boolean (default: ``True``); whether or
+ not to orthonormalize the basis. Doing so is expensive and
+ generally rules out using the rationals as your ``field``, but
+ is required for spectral decompositions.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import random_eja
+
+ TESTS:
+
+ We should compute that an element subalgebra is associative even
+ if we circumvent the element method::
+
+ sage: 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
+
+ @staticmethod
+ def _check_input_field(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")
+
+ @staticmethod
+ def _check_input_axioms(basis, jordan_product, inner_product):
+ 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")
+
+ def __init__(self,
+ basis,
+ jordan_product,
+ inner_product,
+ field=AA,
+ matrix_space=None,
+ orthonormalize=True,
+ associative=None,
+ check_field=True,
+ check_axioms=True,
+ prefix="b"):
+
+ n = len(basis)
+
+ if check_field:
+ self._check_input_field(field)
+
+ 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.
+ self._check_input_axioms(basis, jordan_product, inner_product)
+
+ if n <= 1:
+ # All zero- and one-dimensional algebras are just the real
+ # numbers with (some positive multiples of) the usual
+ # multiplication as its Jordan and inner-product.
+ associative = True
+ if associative is None:
+ # We should figure it out. As with check_axioms, we have to do
+ # this without the help of the _jordan_product_is_associative()
+ # method because we need to know the category before we
+ # initialize the algebra.
+ associative = all( jordan_product(jordan_product(bi,bj),bk)
+ ==
+ jordan_product(bi,jordan_product(bj,bk))
+ for bi in basis
+ for bj in basis
+ for bk in basis)
+
+ category = EuclideanJordanAlgebras(field)
+
+ if associative:
+ # Element subalgebras can take advantage of this.
+ category = category.Associative()
+
+ # 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.
+
+ 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 hold the orthonormal -> unorthonormal
+ # coordinate transformation. Default to an identity matrix of
+ # the appropriate size to avoid special cases for None
+ # everywhere.
+ self._deortho_matrix = matrix.identity(field,n)
+
+ if orthonormalize:
+ # Save a copy of the un-orthonormalized basis for later.
+ # Convert it to ambient V (vector) coordinates while we're
+ # at it, because we'd have to do it later anyway.
+ deortho_vector_basis = tuple( V(_all2list(b)) for b in basis )
+
+ from mjo.eja.eja_utils import gram_schmidt
+ basis = tuple(gram_schmidt(basis, inner_product))
+
+ # Save the (possibly orthonormalized) matrix basis for
+ # later, as well as the space that its elements live in.
+ # In most cases we can deduce the matrix space, but when
+ # n == 0 (that is, there are no basis elements) we cannot.
+ self._matrix_basis = basis
+ if matrix_space is None:
+ self._matrix_space = self._matrix_basis[0].parent()
+ else:
+ self._matrix_space = matrix_space
+
+ # Now create the vector space for the algebra, which will have
+ # its own set of non-ambient coordinates (in terms of the
+ # supplied basis).
+ vector_basis = tuple( V(_all2list(b)) for b in basis )
+
+ # Save the span of our matrix basis (when written out as long
+ # vectors) because otherwise we'll have to reconstruct it
+ # every time we want to coerce a matrix into the algebra.
+ self._matrix_span = V.span_of_basis( vector_basis, check=check_axioms)
+
+ if orthonormalize:
+ # Now "self._matrix_span" is the vector space of our
+ # algebra coordinates. The variables "X0", "X1",... refer
+ # to the entries of vectors in self._matrix_span. Thus to
+ # convert back and forth between the orthonormal
+ # coordinates and the given ones, we need to stick the
+ # original basis in self._matrix_span.
+ U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
+ self._deortho_matrix = matrix.column( 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)
+ zed = self.zero()
+ self._multiplication_table = [ [zed 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 = self._matrix_span.coordinate_vector(V(_all2list(elt)))
+ self._multiplication_table[i][j] = self.from_vector(elt)
+
+ if not orthonormalize:
+ # If we're orthonormalizing the basis with respect
+ # to an inner-product, then the inner-product
+ # matrix with respect to the resulting basis is
+ # just going to be the identity.
+ ip = inner_product(q_i, q_j)
+ self._inner_product_matrix[i,j] = ip
+ self._inner_product_matrix[j,i] = ip
+
+ self._inner_product_matrix._cache = {'hermitian': True}
+ self._inner_product_matrix.set_immutable()
+
+ if check_axioms:
+ if not self._is_jordanian():
+ raise ValueError("Jordan identity does not hold")
+ if not self._inner_product_is_associative():
+ raise ValueError("inner product is not associative")
+
+
def _coerce_map_from_base_ring(self):
"""
Disable the map from the base ring into the algebra.
TESTS::
- sage: set_random_seed()
sage: J = random_eja()
sage: J(1)
Traceback (most recent call last):
"""
return None
- def __init__(self,
- field,
- multiplication_table,
- inner_product_table,
- prefix='e',
- category=None,
- matrix_basis=None,
- check_field=True,
- check_axioms=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: 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
+
"""
- INPUT:
+ # 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]
- * field -- the scalar field for this algebra (must be real)
+ def inner_product(self, x, y):
+ """
+ The inner product associated with this Euclidean Jordan algebra.
- * multiplication_table -- the multiplication table for this
- algebra's implicit basis. Only the lower-triangular portion
- of the table is used, since the multiplication is assumed
- to be commutative.
+ 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: 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: 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
- An error is raised if the Jordan product is not commutative::
+ """
+ B = self._inner_product_matrix
+ return (B*x.to_vector()).inner_product(y.to_vector())
- sage: JP = ((1,2),(0,0))
- sage: IP = ((1,0),(0,1))
- sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
- Traceback (most recent call last):
- ...
- ValueError: Jordan product is not commutative
- An error is raised if the inner-product is not commutative::
+ def is_associative(self):
+ r"""
+ Return whether or not this algebra's Jordan product is associative.
- sage: JP = ((1,0),(0,1))
- sage: IP = ((1,2),(0,0))
- sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
- Traceback (most recent call last):
- ...
- ValueError: inner-product is not commutative
+ 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_field=True``::
+ The values we've presupplied to the constructors agree with
+ the computation::
- sage: JordanSpinEJA(2, field=QQbar)
- Traceback (most recent call last):
- ...
- ValueError: scalar field is not real
- sage: JordanSpinEJA(2, field=QQbar, check_field=False)
- Euclidean Jordan algebra of dimension 2 over Algebraic Field
+ sage: J = random_eja()
+ sage: J.is_associative() == J._jordan_product_is_associative()
+ True
- The multiplication table must be square with ``check_axioms=True``::
+ """
+ R = self.base_ring()
- sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
- 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
- The multiplication and inner-product tables must be the same
- size (and in particular, the inner-product table must also be
- square) with ``check_axioms=True``::
+ 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)
- sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
- Traceback (most recent call last):
- ...
- ValueError: multiplication and inner-product tables are
- different sizes
- sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
- Traceback (most recent call last):
- ...
- ValueError: multiplication and inner-product tables are
- different sizes
+ 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_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")
+ 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
- # The multiplication and inner-product tables should be square
- # if the user wants us to verify them. And we verify them as
- # soon as possible, because we want to exploit their symmetry.
- n = len(multiplication_table)
- if check_axioms:
- if not all( len(l) == n for l in multiplication_table ):
- raise ValueError("multiplication table is not square")
-
- # If the multiplication table is square, we can check if
- # the inner-product table is square by comparing it to the
- # multiplication table's dimensions.
- msg = "multiplication and inner-product tables are different sizes"
- if not len(inner_product_table) == n:
- raise ValueError(msg)
-
- if not all( len(l) == n for l in inner_product_table ):
- raise ValueError(msg)
-
- # Check commutativity of the Jordan product (symmetry of
- # the multiplication table) and the commutativity of the
- # inner-product (symmetry of the inner-product table)
- # first if we're going to check them at all.. This has to
- # be done before we define product_on_basis(), because
- # that method assumes that self._multiplication_table is
- # symmetric. And it has to be done before we build
- # self._inner_product_matrix, because the process used to
- # construct it assumes symmetry as well.
- if not all( multiplication_table[j][i]
- == multiplication_table[i][j]
- for i in range(n)
- for j in range(i+1) ):
- raise ValueError("Jordan product is not commutative")
-
- if not all( inner_product_table[j][i]
- == inner_product_table[i][j]
- for i in range(n)
- for j in range(i+1) ):
- raise ValueError("inner-product is not commutative")
-
- self._matrix_basis = matrix_basis
-
- if category is None:
- category = MagmaticAlgebras(field).FiniteDimensional()
- category = category.WithBasis().Unital()
-
- 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.
- #
- # Note: we take advantage of symmetry here, and only store
- # the lower-triangular portion of the table.
- self._multiplication_table = [ [ self.vector_space().zero()
- for j in range(i+1) ]
- for i in range(n) ]
+ 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)
- for i in range(n):
- for j in range(i+1):
- elt = self.from_vector(multiplication_table[i][j])
- self._multiplication_table[i][j] = elt
-
- self._multiplication_table = tuple(map(tuple, self._multiplication_table))
-
- # Save our inner product as a matrix, since the efficiency of
- # matrix multiplication will usually outweigh the fact that we
- # have to store a redundant upper- or lower-triangular part.
- # Pre-cache the fact that these are Hermitian (real symmetric,
- # in fact) in case some e.g. matrix multiplication routine can
- # take advantage of it.
- ip_matrix_constructor = lambda i,j: inner_product_table[i][j] if j <= i else inner_product_table[j][i]
- self._inner_product_matrix = matrix(field, n, ip_matrix_constructor)
- self._inner_product_matrix._cache = {'hermitian': True}
- self._inner_product_matrix.set_immutable()
+ if diff.abs() > epsilon:
+ return False
- 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")
+ return True
def _element_constructor_(self, elt):
"""
- Construct an element of this algebra from its vector or matrix
- representation.
+ Construct an element of this algebra or a subalgebra from its
+ EJA element, vector, or matrix representation.
This gets called only after the parent element _call_ method
fails to find a coercion for the argument.
SETUP::
- sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: JordanSpinEJA,
....: HadamardEJA,
....: RealSymmetricEJA)
...
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
+
+ Subalgebra elements are embedded into the superalgebra::
+
+ sage: J = JordanSpinEJA(3)
+ sage: J.one()
+ b0
+ sage: x = sum(J.gens())
+ sage: A = x.subalgebra_generated_by()
+ sage: J(A.one())
+ b0
+
TESTS:
- Ensure that we can convert any element of the two non-matrix
- simple algebras (whose matrix representations are columns)
- back and forth faithfully::
+ Ensure that we can convert any element back and forth
+ faithfully between its matrix and algebra representations::
- sage: set_random_seed()
- sage: J = HadamardEJA.random_instance()
- sage: x = J.random_element()
- sage: J(x.to_vector().column()) == x
- True
- sage: J = JordanSpinEJA.random_instance()
+ sage: J = random_eja()
sage: x = J.random_element()
- sage: J(x.to_vector().column()) == x
+ sage: J(x.to_matrix()) == x
True
+ We cannot coerce elements between algebras just because their
+ matrix representations are compatible::
+
+ sage: J1 = HadamardEJA(3)
+ sage: J2 = JordanSpinEJA(3)
+ sage: J2(J1.one())
+ Traceback (most recent call last):
+ ...
+ ValueError: not an element of this algebra
+ sage: J1(J2.zero())
+ Traceback (most recent call last):
+ ...
+ ValueError: not an element of this algebra
+
"""
msg = "not an element of this algebra"
- if elt == 0:
- # The superclass implementation of random_element()
- # needs to be able to coerce "0" into the algebra.
- return self.zero()
- elif elt in self.base_ring():
+ if elt in self.base_ring():
# Ensure that no base ring -> algebra coercion is performed
# by this method. There's some stupidity in sage that would
# otherwise propagate to this method; for example, sage thinks
# that the integer 3 belongs to the space of 2-by-2 matrices.
raise ValueError(msg)
+ if hasattr(elt, 'superalgebra_element'):
+ # Handle subalgebra elements
+ if elt.parent().superalgebra() == self:
+ return elt.superalgebra_element()
+
+ if hasattr(elt, 'sparse_vector'):
+ # Convert a vector into a column-matrix. We check for
+ # "sparse_vector" and not "column" because matrices also
+ # have a "column" method.
+ elt = elt.column()
+
if elt not in self.matrix_space():
raise ValueError(msg)
# closure whereas the base ring of the 3-by-3 identity matrix
# could be QQ instead of QQbar.
#
- # We pass check=False because the matrix basis is "guaranteed"
- # to be linearly independent... right? Ha ha.
- V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
- W = V.span_of_basis( (_mat2vec(s) for s in self.matrix_basis()),
- check=False)
+ # 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.
+ #
+ elt = self._matrix_span.ambient_vector_space()(_all2list(elt))
try:
- coords = W.coordinate_vector(_mat2vec(elt))
+ coords = self._matrix_span.coordinate_vector(elt)
except ArithmeticError: # vector is not in free module
raise ValueError(msg)
fmt = "Euclidean Jordan algebra of dimension {} over {}"
return fmt.format(self.dimension(), self.base_ring())
- def product_on_basis(self, i, j):
- # 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 _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( 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_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 _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 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, field=RDF) with check_axioms=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):
sage: J = JordanSpinEJA(3)
sage: p = J.characteristic_polynomial_of(); p
- X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
+ X0^2 - X1^2 - X2^2 + (-2*t)*X0 + t^2
sage: xvec = J.one().to_vector()
sage: p(*xvec)
t^2 - 2*t + 1
sage: J = HadamardEJA(2)
sage: J.coordinate_polynomial_ring()
- Multivariate Polynomial Ring in X1, X2...
+ Multivariate Polynomial Ring in X0, X1...
sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
sage: J.coordinate_polynomial_ring()
- Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
+ Multivariate Polynomial Ring in X0, X1, X2, X3, X4, X5...
"""
- var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) )
+ var_names = tuple( "X%d" % z for z in range(self.dimension()) )
return PolynomialRing(self.base_ring(), var_names)
def inner_product(self, x, y):
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)
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)
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: 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 |
+----++----+----+----+----+
"""
n = self.dimension()
- M = [ [ self.zero() for j in range(n) ]
- for i in range(n) ]
- for i in range(n):
- for j in range(i+1):
- M[i][j] = self._multiplication_table[i][j]
- M[j][i] = M[i][j]
+ # Prepend the header row.
+ M = [["*"] + list(self.gens())]
- for i in range(n):
- # Prepend the left "header" column entry Can't do this in
- # the loop because it messes up the symmetry.
- M[i] = [self.monomial(i)] + M[i]
+ # 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) ]
- # Prepend the header row.
- M = [["*"] + list(self.gens())] + M
return table(M, header_row=True, header_column=True, frame=True)
Why implement this for non-matrix algebras? Avoiding special
cases for the :class:`BilinearFormEJA` pays with simplicity in
its own right. But mainly, we would like to be able to assume
- that elements of a :class:`DirectSumEJA` can be displayed
+ that elements of a :class:`CartesianProductEJA` can be displayed
nicely, without having to have special classes for direct sums
one of whose components was a matrix algebra.
sage: J = RealSymmetricEJA(2)
sage: J.basis()
- Finite family {0: e0, 1: e1, 2: e2}
+ Finite family {0: b0, 1: b1, 2: b2}
sage: J.matrix_basis()
(
[1 0] [ 0 0.7071067811865475?] [0 0]
sage: J = JordanSpinEJA(2)
sage: J.basis()
- Finite family {0: e0, 1: e1}
+ Finite family {0: b0, 1: b1}
sage: J.matrix_basis()
(
[1] [0]
[0], [1]
)
"""
- if self._matrix_basis is None:
- M = self.matrix_space()
- return tuple( M(b.to_vector()) for b in self.basis() )
- else:
- return self._matrix_basis
+ return self._matrix_basis
def matrix_space(self):
we think of them as matrices (including column vectors of the
appropriate size).
- Generally this will be an `n`-by-`1` column-vector space,
+ "By default" this will be an `n`-by-`1` column-matrix space,
except when the algebra is trivial. There it's `n`-by-`n`
(where `n` is zero), to ensure that two elements of the matrix
- space (empty matrices) can be multiplied.
+ space (empty matrices) can be multiplied. For algebras of
+ matrices, this returns the space in which their
+ real embeddings live.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+ ....: JordanSpinEJA,
+ ....: QuaternionHermitianEJA,
+ ....: TrivialEJA)
+
+ EXAMPLES:
+
+ By default, the matrix representation is just a column-matrix
+ equivalent to the vector representation::
+
+ sage: J = JordanSpinEJA(3)
+ sage: J.matrix_space()
+ Full MatrixSpace of 3 by 1 dense matrices over Algebraic
+ Real Field
+
+ The matrix representation in the trivial algebra is
+ zero-by-zero instead of the usual `n`-by-one::
+
+ sage: J = TrivialEJA()
+ sage: J.matrix_space()
+ Full MatrixSpace of 0 by 0 dense matrices over Algebraic
+ Real Field
+
+ The matrix space for complex/quaternion Hermitian matrix EJA
+ is the space in which their real-embeddings live, not the
+ original complex/quaternion matrix space::
+
+ sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+ sage: J.matrix_space()
+ Module of 2 by 2 matrices with entries in Algebraic Field over
+ the scalar ring Rational Field
+ sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
+ sage: J.matrix_space()
+ Module of 1 by 1 matrices with entries in Quaternion
+ Algebra (-1, -1) with base ring Rational Field over
+ the scalar ring Rational Field
- Matrix algebras override this with something more useful.
"""
- if self.is_trivial():
- return MatrixSpace(self.base_ring(), 0)
- elif self._matrix_basis is None or len(self._matrix_basis) == 0:
- return MatrixSpace(self.base_ring(), self.dimension(), 1)
- else:
- return self._matrix_basis[0].matrix_space()
+ return self._matrix_space
@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
- The matrix of the unit element's operator is the identity::
+ ::
+
+ 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,
+ regardless of the base field and whether or not we
+ orthonormalize::
- sage: set_random_seed()
sage: J = random_eja()
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()
+ sage: actual = A.one().operator().matrix()
+ sage: expected = matrix.identity(A.base_ring(), A.dimension())
+ sage: actual == expected
+ True
+
+ ::
+
+ 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: J = random_eja(field=QQ, orthonormalize=False)
+ sage: cached = J.one()
+ sage: J.one.clear_cache()
+ sage: J.one() == cached
+ True
+
"""
# We can brute-force compute the matrices of the operators
# that correspond to the basis elements of this algebra.
#
# Of course, matrices aren't vectors in sage, so we have to
# appeal to the "long vectors" isometry.
- oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
+
+ V = VectorSpace(self.base_ring(), self.dimension()**2)
+ oper_vecs = [ V(g.operator().matrix().list()) for g in self.gens() ]
# Now we use basic linear algebra to find the coefficients,
# of the matrices-as-vectors-linear-combination, which should
# 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()) )
+ b = V( matrix.identity(self.base_ring(), self.dimension()).list() )
# Now if there's an identity element in the algebra, this
# should work. We solve on the left to avoid having to
Every algebra decomposes trivially with respect to its identity
element::
- sage: set_random_seed()
sage: J = random_eja()
sage: J0,J5,J1 = J.peirce_decomposition(J.one())
sage: J0.dimension() == 0 and J5.dimension() == 0
elements in the two subalgebras are the projections onto their
respective subspaces of the superalgebra's identity element::
- sage: set_random_seed()
sage: J = random_eja()
sage: x = J.random_element()
sage: if not J.is_trivial():
if not c.is_idempotent():
raise ValueError("element is not idempotent: %s" % c)
- from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra
-
# Default these to what they should be if they turn out to be
# trivial, because eigenspaces_left() won't return eigenvalues
# corresponding to trivial spaces (e.g. it returns only the
# eigenspace corresponding to lambda=1 if you take the
# decomposition relative to the identity element).
- trivial = FiniteDimensionalEuclideanJordanSubalgebra(self, ())
+ trivial = self.subalgebra((), check_axioms=False)
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,
- check_axioms=False)
+ subalg = self.subalgebra(gens, check_axioms=False)
if eigval == 0:
J0 = subalg
elif eigval == 1:
for idx in range(count) )
+ def operator_polynomial_matrix(self):
+ r"""
+ Return the matrix of polynomials (over this algebra's
+ :meth:`coordinate_polynomial_ring`) that, when evaluated at
+ the basis coordinates of an element `x`, produces the basis
+ representation of `L_{x}`.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (HadamardEJA,
+ ....: JordanSpinEJA)
+
+ EXAMPLES::
+
+ sage: J = HadamardEJA(4)
+ sage: L_x = J.operator_polynomial_matrix()
+ sage: L_x
+ [X0 0 0 0]
+ [ 0 X1 0 0]
+ [ 0 0 X2 0]
+ [ 0 0 0 X3]
+ sage: x = J.one()
+ sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector())
+ sage: L_x.subs(dict(d))
+ [1 0 0 0]
+ [0 1 0 0]
+ [0 0 1 0]
+ [0 0 0 1]
+
+ ::
+
+ sage: J = JordanSpinEJA(4)
+ sage: L_x = J.operator_polynomial_matrix()
+ sage: L_x
+ [X0 X1 X2 X3]
+ [X1 X0 0 0]
+ [X2 0 X0 0]
+ [X3 0 0 X0]
+ sage: x = J.one()
+ sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector())
+ sage: L_x.subs(dict(d))
+ [1 0 0 0]
+ [0 1 0 0]
+ [0 0 1 0]
+ [0 0 0 1]
+
+ """
+ R = self.coordinate_polynomial_ring()
+
+ def L_x_i_j(i,j):
+ # From a result in my book, these are the entries of the
+ # basis representation of L_x.
+ return sum( v*self.monomial(k).operator().matrix()[i,j]
+ for (k,v) in enumerate(R.gens()) )
+
+ n = self.dimension()
+ return matrix(R, n, n, L_x_i_j)
+
@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: J = random_eja()
+ sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
+ True
+
"""
n = self.dimension()
R = self.coordinate_polynomial_ring()
- vars = R.gens()
F = R.fraction_field()
- def L_x_i_j(i,j):
- # From a result in my book, these are the entries of the
- # basis representation of L_x.
- return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
- for k in range(n) )
-
- L_x = matrix(F, n, n, L_x_i_j)
+ L_x = self.operator_polynomial_matrix()
r = None
if self.rank.is_in_cache():
# 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):
positive integer rank, unless the algebra is trivial in
which case its rank will be zero::
- sage: set_random_seed()
sage: J = random_eja()
sage: r = J.rank()
sage: r in ZZ
Ensure that computing the rank actually works, since the ranks
of all simple algebras are known and will be cached by default::
- sage: set_random_seed() # long time
sage: J = random_eja() # long time
- sage: caches = J.rank() # long time
+ sage: cached = J.rank() # long time
sage: J.rank.clear_cache() # long time
sage: J.rank() == cached # long time
True
return len(self._charpoly_coefficients())
+ def subalgebra(self, basis, **kwargs):
+ r"""
+ Create a subalgebra of this algebra from the given basis.
+ """
+ from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+ return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
+
+
def vector_space(self):
"""
Return the vector space that underlies this algebra.
return self.zero().to_vector().parent().ambient_vector_space()
- Element = FiniteDimensionalEuclideanJordanAlgebraElement
-class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
+class RationalBasisEJA(FiniteDimensionalEJA):
r"""
- New class for algebras whose supplied basis elements have all rational entries.
+ Algebras whose supplied basis elements have all rational entries.
SETUP::
jordan_product,
inner_product,
field=AA,
- orthonormalize=True,
- prefix='e',
- category=None,
check_field=True,
- check_axioms=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 ):
+ # Use _all2list to get the vector coordinates of octonion
+ # entries and not the octonions themselves (which are not
+ # rational).
+ if not all( all(b_i in QQ for b_i in _all2list(b))
+ for b in basis ):
raise TypeError("basis not rational")
- # Temporary(?) hack to ensure that the matrix and vector bases
- # are over the same ring.
- basis = tuple( b.change_ring(field) for b in basis )
-
- n = len(basis)
- vector_basis = basis
-
- from sage.structure.element import is_Matrix
- basis_is_matrices = False
+ super().__init__(basis,
+ jordan_product,
+ inner_product,
+ field=field,
+ check_field=check_field,
+ **kwargs)
- degree = 0
- if n > 0:
- if is_Matrix(basis[0]):
- basis_is_matrices = True
- from mjo.eja.eja_utils import _vec2mat
- vector_basis = tuple( map(_mat2vec,basis) )
- degree = basis[0].nrows()**2
- else:
- degree = basis[0].degree()
-
- V = VectorSpace(field, degree)
-
- # If we were asked to orthonormalize, and if the orthonormal
- # basis is different from the given one, then we also want to
- # compute multiplication and inner-product tables for the
- # deorthonormalized basis. These can be used later to
- # construct a deorthonormalized copy of this algebra over QQ
- # in which several operations are much faster.
self._rational_algebra = None
-
- if orthonormalize:
- if self.base_ring() is not QQ:
- # There's no point in constructing the extra algebra if this
- # one is already rational. If the original basis is rational
- # but normalization would make it irrational, then this whole
- # constructor will just fail anyway as it tries to stick an
- # irrational number into a rational algebra.
- #
- # 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 = RationalBasisEuclideanJordanAlgebra(
- basis,
- jordan_product,
- inner_product,
- field=QQ,
- orthonormalize=False,
- prefix=prefix,
- category=category,
- check_field=False,
- check_axioms=False)
-
- # Compute the deorthonormalized tables before we orthonormalize
- # the given basis. The "check" parameter here guarantees that
- # the basis is linearly-independent.
- W = V.span_of_basis( vector_basis, check=check_axioms)
-
- # 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):
- # given basis w.r.t. ambient coords
- q_i = vector_basis[i]
- q_j = vector_basis[j]
-
- if basis_is_matrices:
- q_i = _vec2mat(q_i)
- q_j = _vec2mat(q_j)
-
- elt = jordan_product(q_i, q_j)
- ip = inner_product(q_i, q_j)
-
- if basis_is_matrices:
- # do another mat2vec because the multiplication
- # table is in terms of vectors
- elt = _mat2vec(elt)
-
- # We overwrite the name "vector_basis" in a second, but never modify it
- # in place, to this effectively makes a copy of it.
- deortho_vector_basis = vector_basis
- self._deortho_matrix = None
-
- if orthonormalize:
- from mjo.eja.eja_utils import gram_schmidt
- if basis_is_matrices:
- vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y))
- vector_basis = gram_schmidt(vector_basis, vector_ip)
- else:
- vector_basis = gram_schmidt(vector_basis, inner_product)
-
- # Normalize the "matrix" basis, too!
- basis = vector_basis
-
- if basis_is_matrices:
- basis = tuple( map(_vec2mat,basis) )
-
- W = V.span_of_basis( vector_basis, check=check_axioms)
-
- # 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 )
-
- # If the superclass constructor is going to verify the
- # symmetry of this table, it has better at least be
- # square...
- if check_axioms:
- mult_table = [ [0 for j in range(n)] for i in range(n) ]
- ip_table = [ [0 for j in range(n)] for i in range(n) ]
- else:
- mult_table = [ [0 for j in range(i+1)] for i in range(n) ]
- ip_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 = vector_basis[i]
- q_j = vector_basis[j]
-
- if basis_is_matrices:
- q_i = _vec2mat(q_i)
- q_j = _vec2mat(q_j)
-
- elt = jordan_product(q_i, q_j)
- ip = inner_product(q_i, q_j)
-
- if basis_is_matrices:
- # do another mat2vec because the multiplication
- # table is in terms of vectors
- elt = _mat2vec(elt)
-
- elt = W.coordinate_vector(elt)
- mult_table[i][j] = elt
- ip_table[i][j] = ip
- if check_axioms:
- # The tables are square if we're verifying that they
- # are commutative.
- mult_table[j][i] = elt
- ip_table[j][i] = ip
-
- if basis_is_matrices:
- for m in basis:
- m.set_immutable()
+ if field is not QQ:
+ # There's no point in constructing the extra algebra if this
+ # one is already rational.
+ #
+ # Note: the same Jordan and inner-products work here,
+ # because they are necessarily defined with respect to
+ # ambient coordinates and not any particular basis.
+ self._rational_algebra = FiniteDimensionalEJA(
+ basis,
+ jordan_product,
+ inner_product,
+ field=QQ,
+ matrix_space=self.matrix_space(),
+ associative=self.is_associative(),
+ orthonormalize=False,
+ check_field=False,
+ check_axioms=False)
+
+ def rational_algebra(self):
+ # Using None as a flag here (rather than just assigning "self"
+ # to self._rational_algebra by default) feels a little bit
+ # more sane to me in a garbage-collected environment.
+ if self._rational_algebra is None:
+ return self
else:
- basis = tuple( x.column() for x in basis )
-
- super().__init__(field,
- mult_table,
- ip_table,
- prefix,
- category,
- basis, # matrix basis
- check_field,
- check_axioms)
+ return self._rational_algebra
@cached_method
def _charpoly_coefficients(self):
sage: J = JordanSpinEJA(3)
sage: J._charpoly_coefficients()
- (X1^2 - X2^2 - X3^2, -2*X1)
+ (X0^2 - X1^2 - X2^2, -2*X0)
sage: a0 = J._charpoly_coefficients()[0]
sage: J.base_ring()
Algebraic Real Field
Algebraic Real Field
"""
- if self.base_ring() is QQ or 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.
- superclass = super(RationalBasisEuclideanJordanAlgebra, self)
- return superclass._charpoly_coefficients()
+ if self.rational_algebra() is self:
+ # Bypass the hijinks if they won't benefit us.
+ 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() )
+ for a_i in self.rational_algebra()._charpoly_coefficients() )
- # Now convert the coordinate variables back to the
+ # Otherwise, convert the coordinate variables back to the
# deorthonormalized ones.
R = self.coordinate_polynomial_ring()
from sage.modules.free_module_element import vector
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 ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
+class ConcreteEJA(FiniteDimensionalEJA):
r"""
A class for the Euclidean Jordan algebras that we know by name.
SETUP::
- sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
+ sage: from mjo.eja.eja_algebra import ConcreteEJA
TESTS:
Our basis is normalized with respect to the algebra's inner
product, unless we specify otherwise::
- sage: set_random_seed()
- sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+ sage: J = ConcreteEJA.random_instance()
sage: all( b.norm() == 1 for b in J.gens() )
True
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 = ConcreteEuclideanJordanAlgebra.random_instance()
+ sage: J = ConcreteEJA.random_instance()
sage: x = J.random_element()
sage: x.operator().is_self_adjoint()
True
"""
@staticmethod
- def _max_random_instance_size():
+ def _max_random_instance_dimension():
+ r"""
+ The maximum dimension of any random instance. Ten dimensions seems
+ to be about the point where everything takes a turn for the
+ worse. And dimension ten (but not nine) allows the 4-by-4 real
+ Hermitian matrices, the 2-by-2 quaternion Hermitian matrices,
+ and the 2-by-2 octonion Hermitian matrices.
+ """
+ return 10
+
+ @staticmethod
+ def _max_random_instance_size(max_dimension):
"""
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 algebra when it is used in a random test case. This size
+ (which can be passed to the algebra's constructor) is itself
+ based on the ``max_dimension`` parameter.
This method must be implemented in each subclass.
"""
raise NotImplementedError
@classmethod
- def random_instance(cls, *args, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
- Return a random instance of this type of algebra.
+ Return a random instance of this type of algebra whose dimension
+ is less than or equal to the lesser of ``max_dimension`` and
+ the value returned by ``_max_random_instance_dimension()``. If
+ the dimension bound is omitted, then only the
+ ``_max_random_instance_dimension()`` is used as a bound.
+
+ This method should be implemented in each subclass.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import ConcreteEJA
+
+ TESTS:
+
+ Both the class bound and the ``max_dimension`` argument are upper
+ bounds on the dimension of the algebra returned::
+
+ sage: from sage.misc.prandom import choice
+ sage: eja_class = choice(ConcreteEJA.__subclasses__())
+ sage: class_max_d = eja_class._max_random_instance_dimension()
+ sage: J = eja_class.random_instance(max_dimension=20,
+ ....: field=QQ,
+ ....: orthonormalize=False)
+ sage: J.dimension() <= class_max_d
+ True
+ sage: J = eja_class.random_instance(max_dimension=2,
+ ....: field=QQ,
+ ....: orthonormalize=False)
+ sage: J.dimension() <= 2
+ True
- This method should be implemented in each subclass.
"""
from sage.misc.prandom import choice
eja_class = choice(cls.__subclasses__())
- # These all bubble up to the RationalBasisEuclideanJordanAlgebra
- # superclass constructor, so any (kw)args valid there are also
- # valid here.
- return eja_class.random_instance(*args, **kwargs)
+ # These all bubble up to the RationalBasisEJA superclass
+ # constructor, so any (kw)args valid there are also valid
+ # here.
+ return eja_class.random_instance(max_dimension, *args, **kwargs)
-class MatrixEuclideanJordanAlgebra:
+class HermitianMatrixEJA(FiniteDimensionalEJA):
@staticmethod
- def real_embed(M):
+ def _denormalized_basis(A):
"""
- Embed the matrix ``M`` into a space of real matrices.
+ Returns a basis for the given Hermitian matrix space.
- The matrix ``M`` can have entries in any field at the moment:
- the real numbers, complex numbers, or quaternions. And although
- they are not a field, we can probably support octonions at some
- point, too. This function returns a real matrix that "acts like"
- the original with respect to matrix multiplication; i.e.
+ Why do we embed these? Basically, because all of numerical linear
+ algebra assumes that you're working with vectors consisting of `n`
+ entries from a field and scalars from the same field. There's no way
+ to tell SageMath that (for example) the vectors contain complex
+ numbers, while the scalar field is real.
- real_embed(M*N) = real_embed(M)*real_embed(N)
+ SETUP::
- """
- raise NotImplementedError
+ sage: from mjo.hurwitz import (ComplexMatrixAlgebra,
+ ....: QuaternionMatrixAlgebra,
+ ....: OctonionMatrixAlgebra)
+ sage: from mjo.eja.eja_algebra import HermitianMatrixEJA
+ TESTS::
+
+ sage: n = ZZ.random_element(1,5)
+ sage: A = MatrixSpace(QQ, n)
+ sage: B = HermitianMatrixEJA._denormalized_basis(A)
+ sage: all( M.is_hermitian() for M in B)
+ True
+
+ ::
+
+ sage: n = ZZ.random_element(1,5)
+ sage: A = ComplexMatrixAlgebra(n, scalars=QQ)
+ sage: B = HermitianMatrixEJA._denormalized_basis(A)
+ sage: all( M.is_hermitian() for M in B)
+ True
+
+ ::
+
+ sage: n = ZZ.random_element(1,5)
+ sage: A = QuaternionMatrixAlgebra(n, scalars=QQ)
+ sage: B = HermitianMatrixEJA._denormalized_basis(A)
+ sage: all( M.is_hermitian() for M in B )
+ True
+
+ ::
+
+ sage: n = ZZ.random_element(1,5)
+ sage: A = OctonionMatrixAlgebra(n, scalars=QQ)
+ sage: B = HermitianMatrixEJA._denormalized_basis(A)
+ sage: all( M.is_hermitian() for M in B )
+ True
- @staticmethod
- def real_unembed(M):
- """
- The inverse of :meth:`real_embed`.
"""
- raise NotImplementedError
+ # These work for real MatrixSpace, whose monomials only have
+ # two coordinates (because the last one would always be "1").
+ es = A.base_ring().gens()
+ gen = lambda A,m: A.monomial(m[:2])
+
+ if hasattr(A, 'entry_algebra_gens'):
+ # We've got a MatrixAlgebra, and its monomials will have
+ # three coordinates.
+ es = A.entry_algebra_gens()
+ gen = lambda A,m: A.monomial(m)
+
+ basis = []
+ for i in range(A.nrows()):
+ for j in range(i+1):
+ if i == j:
+ E_ii = gen(A, (i,j,es[0]))
+ basis.append(E_ii)
+ else:
+ for e in es:
+ E_ij = gen(A, (i,j,e))
+ E_ij += E_ij.conjugate_transpose()
+ basis.append(E_ij)
+
+ return tuple( basis )
@staticmethod
def jordan_product(X,Y):
return (X*Y + Y*X)/2
- @classmethod
- def trace_inner_product(cls,X,Y):
- Xu = cls.real_unembed(X)
- Yu = cls.real_unembed(Y)
- tr = (Xu*Yu).trace()
+ @staticmethod
+ def trace_inner_product(X,Y):
+ r"""
+ A trace inner-product for matrices that aren't embedded in the
+ reals. It takes MATRICES as arguments, not EJA elements.
- try:
- # Works in QQ, AA, RDF, et cetera.
- return tr.real()
- except AttributeError:
- # A quaternion doesn't have a real() method, but does
- # have coefficient_tuple() method that returns the
- # coefficients of 1, i, j, and k -- in that order.
- return tr.coefficient_tuple()[0]
+ SETUP::
+ sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
+ ....: ComplexHermitianEJA,
+ ....: QuaternionHermitianEJA,
+ ....: OctonionHermitianEJA)
-class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
- @staticmethod
- def real_embed(M):
- """
- The identity function, for embedding real matrices into real
- matrices.
- """
- return M
+ EXAMPLES::
+
+ sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False)
+ sage: I = J.one().to_matrix()
+ sage: J.trace_inner_product(I, -I)
+ -2
+
+ ::
+
+ sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+ sage: I = J.one().to_matrix()
+ sage: J.trace_inner_product(I, -I)
+ -2
+
+ ::
+
+ sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
+ sage: I = J.one().to_matrix()
+ sage: J.trace_inner_product(I, -I)
+ -2
+
+ ::
+
+ sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
+ sage: I = J.one().to_matrix()
+ sage: J.trace_inner_product(I, -I)
+ -2
- @staticmethod
- def real_unembed(M):
- """
- The identity function, for unembedding real matrices from real
- matrices.
"""
- return M
+ tr = (X*Y).trace()
+ if hasattr(tr, 'coefficient'):
+ # Works for octonions, and has to come first because they
+ # also have a "real()" method that doesn't return an
+ # element of the scalar ring.
+ return tr.coefficient(0)
+ elif hasattr(tr, 'coefficient_tuple'):
+ # Works for quaternions.
+ return tr.coefficient_tuple()[0]
+
+ # Works for real and complex numbers.
+ return tr.real()
+
+ def __init__(self, matrix_space, **kwargs):
+ # We know this is a valid EJA, but will double-check
+ # if the user passes check_axioms=True.
+ if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
-class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
- RealMatrixEuclideanJordanAlgebra):
+ super().__init__(self._denormalized_basis(matrix_space),
+ self.jordan_product,
+ self.trace_inner_product,
+ field=matrix_space.base_ring(),
+ matrix_space=matrix_space,
+ **kwargs)
+
+ self.rank.set_cache(matrix_space.nrows())
+ self.one.set_cache( self(matrix_space.one()) )
+
+class RealSymmetricEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
"""
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, field=RDF)
+ sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
Euclidean Jordan algebra of dimension 3 over Real Double Field
- sage: RealSymmetricEJA(2, field=RR)
+ sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
Euclidean Jordan algebra of dimension 3 over Real Field with
53 bits of precision
The dimension of this algebra is `(n^2 + n) / 2`::
- sage: set_random_seed()
- sage: n_max = RealSymmetricEJA._max_random_instance_size()
- sage: n = ZZ.random_element(1, n_max)
+ sage: d = RealSymmetricEJA._max_random_instance_dimension()
+ sage: n = RealSymmetricEJA._max_random_instance_size(d)
sage: J = RealSymmetricEJA(n)
sage: J.dimension() == (n^2 + n)/2
True
The Jordan multiplication is what we think it is::
- sage: set_random_seed()
sage: J = RealSymmetricEJA.random_instance()
sage: x,y = J.random_elements(2)
sage: actual = (x*y).to_matrix()
Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
"""
- @classmethod
- def _denormalized_basis(cls, n):
- """
- Return a basis for the space of real symmetric n-by-n matrices.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import RealSymmetricEJA
-
- TESTS::
-
- sage: set_random_seed()
- sage: n = ZZ.random_element(1,5)
- sage: B = RealSymmetricEJA._denormalized_basis(n)
- sage: all( M.is_symmetric() for M in B)
- True
-
- """
- # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
- # coordinates.
- S = []
- for i in range(n):
- for j in range(i+1):
- Eij = matrix(ZZ, n, lambda k,l: k==i and l==j)
- if i == j:
- Sij = Eij
- else:
- Sij = Eij + Eij.transpose()
- S.append(Sij)
- return tuple(S)
-
-
@staticmethod
- def _max_random_instance_size():
- return 4 # Dimension 10
+ def _max_random_instance_size(max_dimension):
+ # Obtained by solving d = (n^2 + n)/2.
+ # The ZZ-int-ZZ thing is just "floor."
+ return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/2 - 1/2))
@classmethod
- def random_instance(cls, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
Return a random instance of this type of algebra.
"""
- n = ZZ.random_element(cls._max_random_instance_size() + 1)
+ class_max_d = cls._max_random_instance_dimension()
+ if (max_dimension is None or max_dimension > class_max_d):
+ max_dimension = class_max_d
+ max_size = cls._max_random_instance_size(max_dimension)
+ n = ZZ.random_element(max_size + 1)
return cls(n, **kwargs)
- def __init__(self, n, **kwargs):
- super(RealSymmetricEJA, self).__init__(self._denormalized_basis(n),
- self.jordan_product,
- self.trace_inner_product,
- **kwargs)
- self.rank.set_cache(n)
- self.one.set_cache(self(matrix.identity(ZZ,n)))
-
-
-class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
- @staticmethod
- def real_embed(M):
- """
- Embed the n-by-n complex matrix ``M`` into the space of real
- matrices of size 2n-by-2n via the map the sends each entry `z = a +
- bi` to the block matrix ``[[a,b],[-b,a]]``.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import \
- ....: ComplexMatrixEuclideanJordanAlgebra
-
- EXAMPLES::
-
- sage: F = QuadraticField(-1, 'I')
- sage: x1 = F(4 - 2*i)
- sage: x2 = F(1 + 2*i)
- sage: x3 = F(-i)
- sage: x4 = F(6)
- sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
- sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
- [ 4 -2| 1 2]
- [ 2 4|-2 1]
- [-----+-----]
- [ 0 -1| 6 0]
- [ 1 0| 0 6]
-
- TESTS:
-
- Embedding is a homomorphism (isomorphism, in fact)::
-
- sage: set_random_seed()
- 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*Ye == XYe
- True
-
- """
- 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]]))
-
- return matrix.block(field, n, blocks)
-
-
- @staticmethod
- def real_unembed(M):
- """
- The inverse of _embed_complex_matrix().
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import \
- ....: ComplexMatrixEuclideanJordanAlgebra
-
- EXAMPLES::
-
- sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
- ....: [-2, 1, -4, 3],
- ....: [ 9, 10, 11, 12],
- ....: [-10, 9, -12, 11] ])
- sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
- [ 2*I + 1 4*I + 3]
- [ 10*I + 9 12*I + 11]
-
- TESTS:
-
- Unembedding is the inverse of embedding::
-
- sage: set_random_seed()
- sage: F = QuadraticField(-1, 'I')
- sage: M = random_matrix(F, 3)
- sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
- sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
- True
-
- """
- n = ZZ(M.nrows())
- if M.ncols() != n:
- raise ValueError("the matrix 'M' must be square")
- if not n.mod(2).is_zero():
- raise ValueError("the matrix 'M' must be a complex embedding")
-
- # If "M" was normalized, its base ring might have roots
- # adjoined and they can stick around after unembedding.
- field = M.base_ring()
- R = PolynomialRing(field, 'z')
- z = R.gen()
- 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())
- 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]
- if submat[0,0] != submat[1,1]:
- raise ValueError('bad on-diagonal submatrix')
- if submat[0,1] != -submat[1,0]:
- raise ValueError('bad off-diagonal submatrix')
- z = submat[0,0] + submat[0,1]*i
- elements.append(z)
-
- return matrix(F, n/2, elements)
-
-
- @classmethod
- def trace_inner_product(cls,X,Y):
- """
- Compute a matrix 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`::
+ def __init__(self, n, field=AA, **kwargs):
+ A = MatrixSpace(field, n)
+ super().__init__(A, **kwargs)
- 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 = ComplexHermitianEJA.real_unembed(Xe)
- sage: Y = ComplexHermitianEJA.real_unembed(Ye)
- sage: expected = (X*Y).trace().real()
- sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye)
- sage: actual == expected
- True
+ from mjo.eja.eja_cache import real_symmetric_eja_coeffs
+ a = real_symmetric_eja_coeffs(self)
+ if a is not None:
+ self.rational_algebra()._charpoly_coefficients.set_cache(a)
- """
- return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2
-class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
- ComplexMatrixEuclideanJordanAlgebra):
+class ComplexHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
"""
The rank-n simple EJA consisting of complex Hermitian n-by-n
matrices over the real numbers, the usual symmetric Jordan product,
EXAMPLES:
- In theory, our "field" can be any subfield of the reals::
+ In theory, our "field" can be any subfield of the reals, but we
+ can't use inexact real fields at the moment because SageMath
+ doesn't know how to convert their elements into complex numbers,
+ or even into algebraic reals::
- sage: ComplexHermitianEJA(2, field=RDF)
- Euclidean Jordan algebra of dimension 4 over Real Double Field
- sage: ComplexHermitianEJA(2, field=RR)
- Euclidean Jordan algebra of dimension 4 over Real Field with
- 53 bits of precision
+ sage: QQbar(RDF(1))
+ Traceback (most recent call last):
+ ...
+ TypeError: Illegal initializer for algebraic number
+ sage: AA(RR(1))
+ Traceback (most recent call last):
+ ...
+ TypeError: Illegal initializer for algebraic number
+
+ This causes the following error when we try to scale a matrix of
+ complex numbers by an inexact real number::
+
+ sage: ComplexHermitianEJA(2,field=RR)
+ Traceback (most recent call last):
+ ...
+ TypeError: Unable to coerce entries (=(1.00000000000000,
+ -0.000000000000000)) to coefficients in Algebraic Real Field
TESTS:
The dimension of this algebra is `n^2`::
- sage: set_random_seed()
- sage: n_max = ComplexHermitianEJA._max_random_instance_size()
- sage: n = ZZ.random_element(1, n_max)
+ sage: d = ComplexHermitianEJA._max_random_instance_dimension()
+ sage: n = ComplexHermitianEJA._max_random_instance_size(d)
sage: J = ComplexHermitianEJA(n)
sage: J.dimension() == n^2
True
The Jordan multiplication is what we think it is::
- sage: set_random_seed()
sage: J = ComplexHermitianEJA.random_instance()
sage: x,y = J.random_elements(2)
sage: actual = (x*y).to_matrix()
Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
"""
+ def __init__(self, n, field=AA, **kwargs):
+ from mjo.hurwitz import ComplexMatrixAlgebra
+ A = ComplexMatrixAlgebra(n, scalars=field)
+ super().__init__(A, **kwargs)
- @classmethod
- def _denormalized_basis(cls, n):
- """
- Returns a basis for the space of complex Hermitian n-by-n matrices.
-
- Why do we embed these? Basically, because all of numerical linear
- algebra assumes that you're working with vectors consisting of `n`
- entries from a field and scalars from the same field. There's no way
- to tell SageMath that (for example) the vectors contain complex
- numbers, while the scalar field is real.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
-
- TESTS::
-
- sage: set_random_seed()
- sage: n = ZZ.random_element(1,5)
- sage: field = QuadraticField(2, 'sqrt2')
- sage: B = ComplexHermitianEJA._denormalized_basis(n)
- sage: all( M.is_symmetric() for M in B)
- True
-
- """
- field = ZZ
- R = PolynomialRing(field, 'z')
- z = R.gen()
- F = field.extension(z**2 + 1, 'I')
- I = F.gen(1)
-
- # This is like the symmetric case, but we need to be careful:
- #
- # * We want conjugate-symmetry, not just symmetry.
- # * The diagonal will (as a result) be real.
- #
- S = []
- for i in range(n):
- for j in range(i+1):
- Eij = matrix(F, n, lambda k,l: k==i and l==j)
- if i == j:
- Sij = cls.real_embed(Eij)
- S.append(Sij)
- else:
- # The second one has a minus because it's conjugated.
- Sij_real = cls.real_embed(Eij + Eij.transpose())
- S.append(Sij_real)
- Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
- S.append(Sij_imag)
-
- # Since we embedded these, we can drop back to the "field" that we
- # started with instead of the complex extension "F".
- return tuple( s.change_ring(field) for s in S )
-
-
- def __init__(self, n, **kwargs):
- super(ComplexHermitianEJA, self).__init__(self._denormalized_basis(n),
- self.jordan_product,
- self.trace_inner_product,
- **kwargs)
- self.rank.set_cache(n)
- # TODO: pre-cache the identity!
+ from mjo.eja.eja_cache import complex_hermitian_eja_coeffs
+ a = complex_hermitian_eja_coeffs(self)
+ if a is not None:
+ self.rational_algebra()._charpoly_coefficients.set_cache(a)
@staticmethod
- def _max_random_instance_size():
- return 3 # Dimension 9
+ def _max_random_instance_size(max_dimension):
+ # Obtained by solving d = n^2.
+ # The ZZ-int-ZZ thing is just "floor."
+ return ZZ(int(ZZ(max_dimension).sqrt()))
@classmethod
- def random_instance(cls, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
Return a random instance of this type of algebra.
"""
- n = ZZ.random_element(cls._max_random_instance_size() + 1)
+ class_max_d = cls._max_random_instance_dimension()
+ if (max_dimension is None or max_dimension > class_max_d):
+ max_dimension = class_max_d
+ max_size = cls._max_random_instance_size(max_dimension)
+ n = ZZ.random_element(max_size + 1)
return cls(n, **kwargs)
-class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
- @staticmethod
- def real_embed(M):
- """
- Embed the n-by-n quaternion matrix ``M`` into the space of real
- matrices of size 4n-by-4n by first sending each quaternion entry `z
- = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
- c+di],[-c + di, a-bi]]`, and then embedding those into a real
- matrix.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import \
- ....: QuaternionMatrixEuclideanJordanAlgebra
-
- EXAMPLES::
-
- sage: Q = QuaternionAlgebra(QQ,-1,-1)
- sage: i,j,k = Q.gens()
- sage: x = 1 + 2*i + 3*j + 4*k
- sage: M = matrix(Q, 1, [[x]])
- sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
- [ 1 2 3 4]
- [-2 1 -4 3]
- [-3 4 1 -2]
- [-4 -3 2 1]
-
- Embedding is a homomorphism (isomorphism, in fact)::
-
- sage: set_random_seed()
- 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*Ye == XYe
- True
-
- """
- 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()
-
- blocks = []
- for z in M.list():
- t = z.coefficient_tuple()
- a = t[0]
- b = t[1]
- c = t[2]
- d = t[3]
- cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
- [-c + d*i, a - b*i]])
- realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM)
- blocks.append(realM)
-
- # We should have real entries by now, so use the realest field
- # we've got for the return value.
- return matrix.block(quaternions.base_ring(), n, blocks)
-
-
- @staticmethod
- def real_unembed(M):
- """
- The inverse of _embed_quaternion_matrix().
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import \
- ....: QuaternionMatrixEuclideanJordanAlgebra
-
- EXAMPLES::
-
- sage: M = matrix(QQ, [[ 1, 2, 3, 4],
- ....: [-2, 1, -4, 3],
- ....: [-3, 4, 1, -2],
- ....: [-4, -3, 2, 1]])
- sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
- [1 + 2*i + 3*j + 4*k]
-
- TESTS:
-
- Unembedding is the inverse of embedding::
-
- sage: set_random_seed()
- sage: Q = QuaternionAlgebra(QQ, -1, -1)
- sage: M = random_matrix(Q, 3)
- sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
- sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
- True
-
- """
- n = ZZ(M.nrows())
- if M.ncols() != n:
- raise ValueError("the matrix 'M' must be square")
- if not n.mod(4).is_zero():
- raise ValueError("the matrix 'M' must be a quaternion embedding")
-
- # Use the base ring of the matrix to ensure that its entries can be
- # multiplied by elements of the quaternion algebra.
- field = M.base_ring()
- Q = QuaternionAlgebra(field,-1,-1)
- i,j,k = Q.gens()
-
- # Go top-left to bottom-right (reading order), converting every
- # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
- # quaternion block.
- elements = []
- for l in 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 trace_inner_product(cls,X,Y):
- """
- Compute a matrix 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.to_matrix()
- sage: Ye = y.to_matrix()
- sage: X = QuaternionHermitianEJA.real_unembed(Xe)
- sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
- sage: expected = (X*Y).trace().coefficient_tuple()[0]
- sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye)
- sage: actual == expected
- True
-
- """
- return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4
-
-
-class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
- QuaternionMatrixEuclideanJordanAlgebra):
+class QuaternionHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
r"""
The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
matrices, the usual symmetric Jordan product, and the
In theory, our "field" can be any subfield of the reals::
- sage: QuaternionHermitianEJA(2, field=RDF)
+ sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
Euclidean Jordan algebra of dimension 6 over Real Double Field
- sage: QuaternionHermitianEJA(2, field=RR)
+ sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
Euclidean Jordan algebra of dimension 6 over Real Field with
53 bits of precision
The dimension of this algebra is `2*n^2 - n`::
- sage: set_random_seed()
- sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
- sage: n = ZZ.random_element(1, n_max)
+ sage: d = QuaternionHermitianEJA._max_random_instance_dimension()
+ sage: n = QuaternionHermitianEJA._max_random_instance_size(d)
sage: J = QuaternionHermitianEJA(n)
sage: J.dimension() == 2*(n^2) - n
True
The Jordan multiplication is what we think it is::
- sage: set_random_seed()
sage: J = QuaternionHermitianEJA.random_instance()
sage: x,y = J.random_elements(2)
sage: actual = (x*y).to_matrix()
Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
"""
- @classmethod
- def _denormalized_basis(cls, n):
+ def __init__(self, n, field=AA, **kwargs):
+ from mjo.hurwitz import QuaternionMatrixAlgebra
+ A = QuaternionMatrixAlgebra(n, scalars=field)
+ super().__init__(A, **kwargs)
+
+ from mjo.eja.eja_cache import quaternion_hermitian_eja_coeffs
+ a = quaternion_hermitian_eja_coeffs(self)
+ if a is not None:
+ self.rational_algebra()._charpoly_coefficients.set_cache(a)
+
+
+
+ @staticmethod
+ def _max_random_instance_size(max_dimension):
+ r"""
+ The maximum rank of a random QuaternionHermitianEJA.
"""
- Returns a basis for the space of quaternion Hermitian n-by-n matrices.
+ # Obtained by solving d = 2n^2 - n.
+ # The ZZ-int-ZZ thing is just "floor."
+ return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4))
- Why do we embed these? Basically, because all of numerical
- linear algebra assumes that you're working with vectors consisting
- of `n` entries from a field and scalars from the same field. There's
- no way to tell SageMath that (for example) the vectors contain
- complex numbers, while the scalar field is real.
+ @classmethod
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
+ """
+ Return a random instance of this type of algebra.
+ """
+ class_max_d = cls._max_random_instance_dimension()
+ if (max_dimension is None or max_dimension > class_max_d):
+ max_dimension = class_max_d
+ max_size = cls._max_random_instance_size(max_dimension)
+ n = ZZ.random_element(max_size + 1)
+ return cls(n, **kwargs)
- SETUP::
+class OctonionHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
+ r"""
+ SETUP::
- sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
+ sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
+ ....: OctonionHermitianEJA)
+ sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra
- TESTS::
+ EXAMPLES:
- sage: set_random_seed()
- sage: n = ZZ.random_element(1,5)
- sage: B = QuaternionHermitianEJA._denormalized_basis(n)
- sage: all( M.is_symmetric() for M in B )
- True
+ The 3-by-3 algebra satisfies the axioms of an EJA::
+
+ sage: OctonionHermitianEJA(3, # long time
+ ....: field=QQ, # long time
+ ....: orthonormalize=False, # long time
+ ....: check_axioms=True) # long time
+ Euclidean Jordan algebra of dimension 27 over Rational Field
+
+ After a change-of-basis, the 2-by-2 algebra has the same
+ multiplication table as the ten-dimensional Jordan spin algebra::
+
+ sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ)
+ sage: b = OctonionHermitianEJA._denormalized_basis(A)
+ sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
+ sage: jp = OctonionHermitianEJA.jordan_product
+ sage: ip = OctonionHermitianEJA.trace_inner_product
+ sage: J = FiniteDimensionalEJA(basis,
+ ....: jp,
+ ....: ip,
+ ....: field=QQ,
+ ....: orthonormalize=False)
+ sage: J.multiplication_table()
+ +----++----+----+----+----+----+----+----+----+----+----+
+ | * || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
+ +====++====+====+====+====+====+====+====+====+====+====+
+ | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
+ +----++----+----+----+----+----+----+----+----+----+----+
+ | b1 || b1 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
+ +----++----+----+----+----+----+----+----+----+----+----+
+ | b2 || b2 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
+ +----++----+----+----+----+----+----+----+----+----+----+
+ | b3 || b3 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 |
+ +----++----+----+----+----+----+----+----+----+----+----+
+ | b4 || b4 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 |
+ +----++----+----+----+----+----+----+----+----+----+----+
+ | b5 || b5 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 |
+ +----++----+----+----+----+----+----+----+----+----+----+
+ | b6 || b6 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 |
+ +----++----+----+----+----+----+----+----+----+----+----+
+ | b7 || b7 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 |
+ +----++----+----+----+----+----+----+----+----+----+----+
+ | b8 || b8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 |
+ +----++----+----+----+----+----+----+----+----+----+----+
+ | b9 || b9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 |
+ +----++----+----+----+----+----+----+----+----+----+----+
- """
- field = ZZ
- Q = QuaternionAlgebra(QQ,-1,-1)
- I,J,K = Q.gens()
+ TESTS:
- # This is like the symmetric case, but we need to be careful:
- #
- # * We want conjugate-symmetry, not just symmetry.
- # * The diagonal will (as a result) be real.
- #
- S = []
- for i in range(n):
- for j in range(i+1):
- Eij = matrix(Q, n, lambda k,l: k==i and l==j)
- if i == j:
- Sij = cls.real_embed(Eij)
- S.append(Sij)
- else:
- # The second, third, and fourth ones have a minus
- # because they're conjugated.
- Sij_real = cls.real_embed(Eij + Eij.transpose())
- S.append(Sij_real)
- Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
- S.append(Sij_I)
- Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
- S.append(Sij_J)
- Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
- S.append(Sij_K)
-
- # Since we embedded these, we can drop back to the "field" that we
- # started with instead of the quaternion algebra "Q".
- return tuple( s.change_ring(field) for s in S )
-
-
- def __init__(self, n, **kwargs):
- super(QuaternionHermitianEJA, self).__init__(self._denormalized_basis(n),
- self.jordan_product,
- self.trace_inner_product,
- **kwargs)
- self.rank.set_cache(n)
- # TODO: cache one()!
+ We can actually construct the 27-dimensional Albert algebra,
+ and we get the right unit element if we recompute it::
+
+ sage: J = OctonionHermitianEJA(3, # long time
+ ....: field=QQ, # long time
+ ....: orthonormalize=False) # long time
+ sage: J.one.clear_cache() # long time
+ sage: J.one() # long time
+ b0 + b9 + b26
+ sage: J.one().to_matrix() # long time
+ +----+----+----+
+ | e0 | 0 | 0 |
+ +----+----+----+
+ | 0 | e0 | 0 |
+ +----+----+----+
+ | 0 | 0 | e0 |
+ +----+----+----+
+
+ The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
+ spin algebra, but just to be sure, we recompute its rank::
+
+ sage: J = OctonionHermitianEJA(2, # long time
+ ....: field=QQ, # long time
+ ....: orthonormalize=False) # long time
+ sage: J.rank.clear_cache() # long time
+ sage: J.rank() # long time
+ 2
+ """
@staticmethod
- def _max_random_instance_size():
+ def _max_random_instance_size(max_dimension):
r"""
- The maximum rank of a random QuaternionHermitianEJA.
- """
- return 2 # Dimension 6
+ The maximum rank of a random OctonionHermitianEJA.
+ """
+ # There's certainly a formula for this, but with only four
+ # cases to worry about, I'm not that motivated to derive it.
+ if max_dimension >= 27:
+ return 3
+ elif max_dimension >= 10:
+ return 2
+ elif max_dimension >= 1:
+ return 1
+ else:
+ return 0
@classmethod
- def random_instance(cls, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
Return a random instance of this type of algebra.
"""
- n = ZZ.random_element(cls._max_random_instance_size() + 1)
+ class_max_d = cls._max_random_instance_dimension()
+ if (max_dimension is None or max_dimension > class_max_d):
+ max_dimension = class_max_d
+ max_size = cls._max_random_instance_size(max_dimension)
+ n = ZZ.random_element(max_size + 1)
return cls(n, **kwargs)
+ def __init__(self, n, field=AA, **kwargs):
+ if n > 3:
+ # Otherwise we don't get an EJA.
+ raise ValueError("n cannot exceed 3")
+
+ # We know this is a valid EJA, but will double-check
+ # if the user passes check_axioms=True.
+ if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+ from mjo.hurwitz import OctonionMatrixAlgebra
+ A = OctonionMatrixAlgebra(n, scalars=field)
+ super().__init__(A, **kwargs)
+
+ from mjo.eja.eja_cache import octonion_hermitian_eja_coeffs
+ a = octonion_hermitian_eja_coeffs(self)
+ if a is not None:
+ self.rational_algebra()._charpoly_coefficients.set_cache(a)
+
+
+class AlbertEJA(OctonionHermitianEJA):
+ r"""
+ The Albert algebra is the algebra of three-by-three Hermitian
+ matrices whose entries are octonions.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import AlbertEJA
+
+ EXAMPLES::
+
+ sage: AlbertEJA(field=QQ, orthonormalize=False)
+ Euclidean Jordan algebra of dimension 27 over Rational Field
+ sage: AlbertEJA() # long time
+ Euclidean Jordan algebra of dimension 27 over Algebraic Real Field
+
+ """
+ def __init__(self, *args, **kwargs):
+ super().__init__(3, *args, **kwargs)
+
-class HadamardEJA(ConcreteEuclideanJordanAlgebra):
+class HadamardEJA(RationalBasisEJA, ConcreteEJA):
"""
- Return the Euclidean Jordan Algebra corresponding to the set
- `R^n` under the Hadamard product.
+ Return the Euclidean Jordan algebra on `R^n` with the Hadamard
+ (pointwise real-number multiplication) Jordan product and the
+ usual inner-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.
+ This is nothing more than the Cartesian product of ``n`` copies of
+ the one-dimensional Jordan spin algebra, and is the most common
+ example of a non-simple Euclidean Jordan algebra.
SETUP::
This multiplication table can be verified by hand::
sage: J = HadamardEJA(3)
- sage: e0,e1,e2 = J.gens()
- sage: e0*e0
- e0
- sage: e0*e1
+ sage: b0,b1,b2 = J.gens()
+ sage: b0*b0
+ b0
+ sage: b0*b1
0
- sage: e0*e2
+ sage: b0*b2
0
- sage: e1*e1
- e1
- sage: e1*e2
+ sage: b1*b1
+ b1
+ sage: b1*b2
0
- sage: e2*e2
- e2
+ sage: b2*b2
+ b2
TESTS:
sage: HadamardEJA(3, prefix='r').gens()
(r0, r1, r2)
-
"""
- def __init__(self, n, **kwargs):
- def jordan_product(x,y):
- P = x.parent()
- return P(tuple( xi*yi for (xi,yi) in zip(x,y) ))
- def inner_product(x,y):
- return x.inner_product(y)
-
- # Don't orthonormalize because our basis is already
- # orthonormal with respect to our inner-product.
- if not 'orthonormalize' in kwargs:
- kwargs['orthonormalize'] = False
-
- # But also don't pass check_field=False here, because the user
- # can pass in a field!
- standard_basis = FreeModule(ZZ, n).basis()
- super(HadamardEJA, self).__init__(standard_basis,
- jordan_product,
- inner_product,
- check_axioms=False,
- **kwargs)
- self.rank.set_cache(n)
+ def __init__(self, n, field=AA, **kwargs):
+ MS = MatrixSpace(field, n, 1)
if n == 0:
- self.one.set_cache( self.zero() )
+ jordan_product = lambda x,y: x
+ inner_product = lambda x,y: x
else:
- self.one.set_cache( sum(self.gens()) )
+ def jordan_product(x,y):
+ return MS( xi*yi for (xi,yi) in zip(x,y) )
+
+ def inner_product(x,y):
+ return (x.T*y)[0,0]
+
+ # New defaults for keyword arguments. Don't orthonormalize
+ # because our basis is already orthonormal with respect to our
+ # inner-product. Don't check the axioms, because we know this
+ # is a valid EJA... but do double-check if the user passes
+ # check_axioms=True. Note: we DON'T override the "check_field"
+ # default here, because the user can pass in a field!
+ if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
+ if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+ column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() )
+ super().__init__(column_basis,
+ jordan_product,
+ inner_product,
+ field=field,
+ matrix_space=MS,
+ associative=True,
+ **kwargs)
+ self.rank.set_cache(n)
+
+ self.one.set_cache( self.sum(self.gens()) )
@staticmethod
- def _max_random_instance_size():
+ def _max_random_instance_dimension():
r"""
- The maximum dimension of a random HadamardEJA.
+ There's no reason to go higher than five here. That's
+ enough to get the point across.
"""
return 5
+ @staticmethod
+ def _max_random_instance_size(max_dimension):
+ r"""
+ The maximum size (=dimension) of a random HadamardEJA.
+ """
+ return max_dimension
+
@classmethod
- def random_instance(cls, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
Return a random instance of this type of algebra.
"""
- n = ZZ.random_element(cls._max_random_instance_size() + 1)
+ class_max_d = cls._max_random_instance_dimension()
+ if (max_dimension is None or max_dimension > class_max_d):
+ max_dimension = class_max_d
+ max_size = cls._max_random_instance_size(max_dimension)
+ n = ZZ.random_element(max_size + 1)
return cls(n, **kwargs)
-class BilinearFormEJA(ConcreteEuclideanJordanAlgebra):
+class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):
r"""
The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
with the half-trace inner product and jordan product ``x*y =
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: B11 = matrix.identity(QQ,1)
....: for j in range(n-1) ]
sage: actual == expected
True
+
"""
- def __init__(self, B, **kwargs):
- if not B.is_positive_definite():
- raise ValueError("bilinear form is not positive-definite")
+ def __init__(self, B, field=AA, **kwargs):
+ # The matrix "B" is supplied by the user in most cases,
+ # so it makes sense to check whether or not its positive-
+ # definite unless we are specifically asked not to...
+ if ("check_axioms" not in kwargs) or kwargs["check_axioms"]:
+ if not B.is_positive_definite():
+ raise ValueError("bilinear form is not positive-definite")
+
+ # However, all of the other data for this EJA is computed
+ # by us in manner that guarantees the axioms are
+ # satisfied. So, again, unless we are specifically asked to
+ # verify things, we'll skip the rest of the checks.
+ if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+ n = B.nrows()
+ MS = MatrixSpace(field, n, 1)
def inner_product(x,y):
- return (B*x).inner_product(y)
+ return (y.T*B*x)[0,0]
def jordan_product(x,y):
- P = x.parent()
- x0 = x[0]
- xbar = x[1:]
- y0 = y[0]
- ybar = y[1:]
- z0 = inner_product(x,y)
+ x0 = x[0,0]
+ xbar = x[1:,0]
+ y0 = y[0,0]
+ ybar = y[1:,0]
+ z0 = inner_product(y,x)
zbar = y0*xbar + x0*ybar
- return P((z0,) + tuple(zbar))
+ return MS([z0] + zbar.list())
- n = B.nrows()
- standard_basis = FreeModule(ZZ, n).basis()
- super(BilinearFormEJA, self).__init__(standard_basis,
- jordan_product,
- inner_product,
- **kwargs)
+ column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() )
+
+ # TODO: I haven't actually checked this, but it seems legit.
+ associative = False
+ if n <= 2:
+ associative = True
+
+ super().__init__(column_basis,
+ jordan_product,
+ inner_product,
+ field=field,
+ matrix_space=MS,
+ associative=associative,
+ **kwargs)
# The rank of this algebra is two, unless we're in a
# one-dimensional ambient space (because the rank is bounded
# by the ambient dimension).
self.rank.set_cache(min(n,2))
-
if n == 0:
self.one.set_cache( self.zero() )
else:
self.one.set_cache( self.monomial(0) )
@staticmethod
- def _max_random_instance_size():
+ def _max_random_instance_dimension():
r"""
- The maximum dimension of a random BilinearFormEJA.
+ There's no reason to go higher than five here. That's
+ enough to get the point across.
"""
return 5
+ @staticmethod
+ def _max_random_instance_size(max_dimension):
+ r"""
+ The maximum size (=dimension) of a random BilinearFormEJA.
+ """
+ return max_dimension
+
@classmethod
- def random_instance(cls, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
Return a random instance of this algebra.
"""
- n = ZZ.random_element(cls._max_random_instance_size() + 1)
+ class_max_d = cls._max_random_instance_dimension()
+ if (max_dimension is None or max_dimension > class_max_d):
+ max_dimension = class_max_d
+ max_size = cls._max_random_instance_size(max_dimension)
+ n = ZZ.random_element(max_size + 1)
+
if n.is_zero():
B = matrix.identity(ZZ, n)
return cls(B, **kwargs)
alpha = ZZ.zero()
while alpha.is_zero():
alpha = ZZ.random_element().abs()
+
B22 = M.transpose()*M + alpha*I
from sage.matrix.special import block_matrix
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::
Ensure that we have the usual inner product on `R^n`::
- sage: set_random_seed()
sage: J = JordanSpinEJA.random_instance()
sage: x,y = J.random_elements(2)
sage: actual = x.inner_product(y)
True
"""
- def __init__(self, n, **kwargs):
+ def __init__(self, n, *args, **kwargs):
# This is a special case of the BilinearFormEJA with the
# identity matrix as its bilinear form.
B = matrix.identity(ZZ, n)
# Don't orthonormalize because our basis is already
# orthonormal with respect to our inner-product.
- if not 'orthonormalize' in kwargs:
- kwargs['orthonormalize'] = False
+ 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(JordanSpinEJA, self).__init__(B,
- check_axioms=False,
- **kwargs)
-
- @staticmethod
- def _max_random_instance_size():
- r"""
- The maximum dimension of a random JordanSpinEJA.
- """
- return 5
+ super().__init__(B, *args, **kwargs)
@classmethod
- def random_instance(cls, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **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)
+ class_max_d = cls._max_random_instance_dimension()
+ if (max_dimension is None or max_dimension > class_max_d):
+ max_dimension = class_max_d
+ max_size = cls._max_random_instance_size(max_dimension)
+ n = ZZ.random_element(max_size + 1)
return cls(n, **kwargs)
-class TrivialEJA(ConcreteEuclideanJordanAlgebra):
+class TrivialEJA(RationalBasisEJA, ConcreteEJA):
"""
The trivial Euclidean Jordan algebra consisting of only a zero element.
0
"""
- def __init__(self, **kwargs):
+ def __init__(self, field=AA, **kwargs):
jordan_product = lambda x,y: x
- inner_product = lambda x,y: 0
+ inner_product = lambda x,y: field.zero()
basis = ()
- super(TrivialEJA, self).__init__(basis,
- jordan_product,
- inner_product,
- **kwargs)
+ MS = MatrixSpace(field,0)
+
+ # New defaults for keyword arguments
+ if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
+ if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+ super().__init__(basis,
+ jordan_product,
+ inner_product,
+ associative=True,
+ field=field,
+ matrix_space=MS,
+ **kwargs)
+
# The rank is zero using my definition, namely the dimension of the
# largest subalgebra generated by any element.
self.rank.set_cache(0)
self.one.set_cache( self.zero() )
@classmethod
- def random_instance(cls, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
# We don't take a "size" argument so the superclass method is
- # inappropriate for us.
+ # inappropriate for us. The ``max_dimension`` argument is
+ # included so that if this method is called generically with a
+ # ``max_dimension=<whatever>`` argument, we don't try to pass
+ # it on to the algebra constructor.
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 (random_eja,
+ ....: CartesianProductEJA,
+ ....: ComplexHermitianEJA,
....: HadamardEJA,
- ....: RealSymmetricEJA,
- ....: DirectSumEJA)
+ ....: JordanSpinEJA,
+ ....: RealSymmetricEJA)
- EXAMPLES::
+ EXAMPLES:
+
+ The Jordan product is inherited from our factors and implemented by
+ our CombinatorialFreeModule Cartesian product superclass::
sage: J1 = HadamardEJA(2)
- sage: J2 = RealSymmetricEJA(3)
- sage: J = DirectSumEJA(J1,J2)
- sage: J.dimension()
- 8
+ sage: J2 = RealSymmetricEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: x,y = J.random_elements(2)
+ sage: x*y in J
+ True
+
+ The ability to retrieve the original factors is implemented by our
+ CombinatorialFreeModule Cartesian product superclass::
+
+ sage: J1 = HadamardEJA(2, field=QQ)
+ sage: J2 = JordanSpinEJA(3, field=QQ)
+ sage: J = cartesian_product([J1,J2])
+ sage: J.cartesian_factors()
+ (Euclidean Jordan algebra of dimension 2 over Rational Field,
+ Euclidean Jordan algebra of dimension 3 over Rational Field)
+
+ You can provide more than two factors::
+
+ sage: J1 = HadamardEJA(2)
+ sage: J2 = JordanSpinEJA(3)
+ sage: J3 = RealSymmetricEJA(3)
+ sage: cartesian_product([J1,J2,J3])
+ Euclidean Jordan algebra of dimension 2 over Algebraic Real
+ Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
+ Real Field (+) Euclidean Jordan algebra of dimension 6 over
+ Algebraic Real Field
+
+ Rank is additive on a Cartesian product::
+
+ sage: J1 = HadamardEJA(1)
+ sage: J2 = RealSymmetricEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: J1.rank.clear_cache()
+ sage: J2.rank.clear_cache()
+ sage: J.rank.clear_cache()
sage: J.rank()
- 5
+ 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 |
+ +----++----+----+----+
+
+ The "matrix space" of a Cartesian product always consists of
+ ordered pairs (or triples, or...) whose components are the
+ matrix spaces of its factors::
+
+ sage: J1 = HadamardEJA(2)
+ sage: J2 = ComplexHermitianEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: J.matrix_space()
+ The Cartesian product of (Full MatrixSpace of 2 by 1 dense
+ matrices over Algebraic Real Field, Module of 2 by 2 matrices
+ with entries in Algebraic Field over the scalar ring Algebraic
+ Real Field)
+ sage: J.one().to_matrix()[0]
+ [1]
+ [1]
+ sage: J.one().to_matrix()[1]
+ +---+---+
+ | 1 | 0 |
+ +---+---+
+ | 0 | 1 |
+ +---+---+
TESTS:
- The external direct sum construction is only valid when the two factors
- have the same base ring; an error is raised otherwise::
+ All factors must share the same base field::
- sage: set_random_seed()
- sage: J1 = random_eja(field=AA)
- sage: J2 = random_eja(field=QQ,orthonormalize=False)
- sage: J = DirectSumEJA(J1,J2)
+ sage: J1 = HadamardEJA(2, field=QQ)
+ sage: J2 = RealSymmetricEJA(2)
+ sage: CartesianProductEJA((J1,J2))
Traceback (most recent call last):
...
- ValueError: algebras must share the same base field
+ ValueError: all factors must share the same base field
+ The cached unit element is the same one that would be computed::
+
+ sage: J1 = random_eja() # long time
+ sage: J2 = random_eja() # long time
+ sage: J = cartesian_product([J1,J2]) # long time
+ sage: actual = J.one() # long time
+ sage: J.one.clear_cache() # long time
+ sage: expected = J.one() # long time
+ sage: actual == expected # long time
+ True
"""
- def __init__(self, J1, J2, **kwargs):
- if J1.base_ring() != J2.base_ring():
- raise ValueError("algebras must share the same base field")
- field = J1.base_ring()
-
- self._factors = (J1, J2)
- n1 = J1.dimension()
- n2 = J2.dimension()
- n = n1+n2
- V = VectorSpace(field, n)
- mult_table = [ [ V.zero() for j in range(i+1) ]
- for i in range(n) ]
- for i in range(n1):
- for j in range(i+1):
- p = (J1.monomial(i)*J1.monomial(j)).to_vector()
- mult_table[i][j] = V(p.list() + [field.zero()]*n2)
+ Element = CartesianProductEJAElement
+ 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")
+
+ # Figure out the category to use.
+ associative = all( f.is_associative() for f in factors )
+ category = EuclideanJordanAlgebras(field)
+ if associative: category = category.Associative()
+ category = category.join([category, category.CartesianProducts()])
+
+ # Compute my matrix space. We don't simply use the
+ # ``cartesian_product()`` functor here because it acts
+ # differently on SageMath MatrixSpaces and our custom
+ # MatrixAlgebras, which are CombinatorialFreeModules. We
+ # always want the result to be represented (and indexed) as an
+ # ordered tuple. This category isn't perfect, but is good
+ # enough for what we need to do.
+ MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
+ MS_cat = MS_cat.Unital().CartesianProducts()
+ MS_factors = tuple( J.matrix_space() for J in factors )
+ from sage.sets.cartesian_product import CartesianProduct
+ self._matrix_space = CartesianProduct(MS_factors, MS_cat)
+
+ self._matrix_basis = []
+ zero = self._matrix_space.zero()
+ for i in range(m):
+ for b in factors[i].matrix_basis():
+ z = list(zero)
+ z[i] = b
+ self._matrix_basis.append(z)
+
+ self._matrix_basis = tuple( self._matrix_space(b)
+ for b in self._matrix_basis )
+ n = len(self._matrix_basis)
+
+ # We already have what we need for the super-superclass constructor.
+ CombinatorialFreeModule.__init__(self,
+ field,
+ range(n),
+ prefix="b",
+ category=category,
+ bracket=False)
+
+ # Now create the vector space for the algebra, which will have
+ # its own set of non-ambient coordinates (in terms of the
+ # supplied basis).
+ degree = sum( f._matrix_span.ambient_vector_space().degree()
+ for f in factors )
+ V = VectorSpace(field, degree)
+ vector_basis = tuple( V(_all2list(b)) for b in self._matrix_basis )
+
+ # Save the span of our matrix basis (when written out as long
+ # vectors) because otherwise we'll have to reconstruct it
+ # every time we want to coerce a matrix into the algebra.
+ self._matrix_span = V.span_of_basis( vector_basis, check=False)
+
+ # Since we don't (re)orthonormalize the basis, the FDEJA
+ # constructor is going to set self._deortho_matrix to the
+ # identity matrix. Here we set it to the correct value using
+ # the deortho matrices from our factors.
+ self._deortho_matrix = matrix.block_diagonal(
+ [J._deortho_matrix for J in factors]
+ )
- for i in range(n2):
- for j in range(i+1):
- p = (J2.monomial(i)*J2.monomial(j)).to_vector()
- mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
-
- # TODO: build the IP table here from the two constituent IP
- # matrices (it'll be block diagonal, I think).
- ip_table = [ [ field.zero() for j in range(i+1) ]
- for i in range(n) ]
- super(DirectSumEJA, self).__init__(field,
- mult_table,
- ip_table,
- check_axioms=False,
- **kwargs)
- self.rank.set_cache(J1.rank() + J2.rank())
-
-
- def factors(self):
+ self._inner_product_matrix = matrix.block_diagonal(
+ [J._inner_product_matrix for J in factors]
+ )
+ self._inner_product_matrix._cache = {'hermitian': True}
+ self._inner_product_matrix.set_immutable()
+
+ # Building the multiplication table is a bit more tricky
+ # because we have to embed the entries of the factors'
+ # multiplication tables into the product EJA.
+ zed = self.zero()
+ self._multiplication_table = [ [zed for j in range(i+1)]
+ for i in range(n) ]
+
+ # Keep track of an offset that tallies the dimensions of all
+ # previous factors. If the second factor is dim=2 and if the
+ # first one is dim=3, then we want to skip the first 3x3 block
+ # when copying the multiplication table for the second factor.
+ offset = 0
+ for f in range(m):
+ phi_f = self.cartesian_embedding(f)
+ factor_dim = factors[f].dimension()
+ for i in range(factor_dim):
+ for j in range(i+1):
+ f_ij = factors[f]._multiplication_table[i][j]
+ e = phi_f(f_ij)
+ self._multiplication_table[offset+i][offset+j] = e
+ offset += factor_dim
+
+ self.rank.set_cache(sum(J.rank() for J in factors))
+ ones = tuple(J.one().to_matrix() for J in factors)
+ self.one.set_cache(self(ones))
+
+ def _sets_keys(self):
r"""
- Return the pair of this algebra's factors.
SETUP::
- sage: from mjo.eja.eja_algebra import (HadamardEJA,
- ....: JordanSpinEJA,
- ....: DirectSumEJA)
+ sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+ ....: RealSymmetricEJA)
- EXAMPLES::
+ TESTS:
+
+ The superclass uses ``_sets_keys()`` to implement its
+ ``cartesian_factors()`` method::
- sage: J1 = HadamardEJA(2, field=QQ)
- sage: J2 = JordanSpinEJA(3, field=QQ)
- sage: J = DirectSumEJA(J1,J2)
- sage: J.factors()
- (Euclidean Jordan algebra of dimension 2 over Rational Field,
- Euclidean Jordan algebra of dimension 3 over Rational Field)
+ sage: J1 = RealSymmetricEJA(2,
+ ....: field=QQ,
+ ....: orthonormalize=False,
+ ....: prefix="a")
+ sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+ sage: J = cartesian_product([J1,J2])
+ sage: x = sum(i*J.gens()[i] for i in range(len(J.gens())))
+ sage: x.cartesian_factors()
+ (a1 + 2*a2, 3*b0 + 4*b1 + 5*b2 + 6*b3)
"""
- return self._factors
+ # Copy/pasted from CombinatorialFreeModule_CartesianProduct,
+ # but returning a tuple instead of a list.
+ return tuple(range(len(self.cartesian_factors())))
- def projections(self):
+ def cartesian_factors(self):
+ # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+ return self._sets
+
+ def cartesian_factor(self, i):
r"""
- Return a pair of projections onto this algebra's factors.
+ 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)
+
+ @cached_method
+ def cartesian_projection(self, i):
+ r"""
SETUP::
- sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
- ....: ComplexHermitianEJA,
- ....: DirectSumEJA)
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: JordanSpinEJA,
+ ....: HadamardEJA,
+ ....: RealSymmetricEJA,
+ ....: ComplexHermitianEJA)
- EXAMPLES::
+ 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 = DirectSumEJA(J1,J2)
- sage: (pi_left, pi_right) = J.projections()
- sage: J.one().to_vector()
- (1, 0, 1, 0, 0, 1)
+ 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: 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
"""
- (J1,J2) = self.factors()
- m = J1.dimension()
- n = J2.dimension()
- V_basis = self.vector_space().basis()
- # Need to specify the dimensions explicitly so that we don't
- # wind up with a zero-by-zero matrix when we want e.g. a
- # zero-by-two matrix (important for composing things).
- P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
- P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
- pi_left = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J1,P1)
- pi_right = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J2,P2)
- return (pi_left, pi_right)
-
- def inclusions(self):
- r"""
- Return the pair of inclusion maps from our factors into us.
+ 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,
- ....: RealSymmetricEJA,
- ....: DirectSumEJA)
+ ....: HadamardEJA,
+ ....: RealSymmetricEJA)
- EXAMPLES::
+ EXAMPLES:
+
+ The embedding morphisms are Euclidean Jordan algebra
+ operators::
+
+ sage: J1 = HadamardEJA(2)
+ sage: J2 = RealSymmetricEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: J.cartesian_embedding(0)
+ Linear operator between finite-dimensional Euclidean Jordan
+ algebras represented by the matrix:
+ [1 0]
+ [0 1]
+ [0 0]
+ [0 0]
+ [0 0]
+ Domain: Euclidean Jordan algebra of dimension 2 over
+ Algebraic Real Field
+ Codomain: Euclidean Jordan algebra of dimension 2 over
+ Algebraic Real Field (+) Euclidean Jordan algebra of
+ dimension 3 over Algebraic Real Field
+ sage: J.cartesian_embedding(1)
+ Linear operator between finite-dimensional Euclidean Jordan
+ algebras represented by the matrix:
+ [0 0 0]
+ [0 0 0]
+ [1 0 0]
+ [0 1 0]
+ [0 0 1]
+ Domain: Euclidean Jordan algebra of dimension 3 over
+ Algebraic Real Field
+ Codomain: Euclidean Jordan algebra of dimension 2 over
+ Algebraic Real Field (+) Euclidean Jordan algebra of
+ dimension 3 over Algebraic Real Field
+
+ The embeddings work the way you'd expect on the vector
+ representation of an element::
sage: J1 = JordanSpinEJA(3)
sage: J2 = RealSymmetricEJA(2)
- sage: J = DirectSumEJA(J1,J2)
- sage: (iota_left, iota_right) = J.inclusions()
+ 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()
TESTS:
+ The answer never changes::
+
+ 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 = DirectSumEJA(J1,J2)
- sage: (iota_left, iota_right) = J.inclusions()
- sage: (pi_left, pi_right) = J.projections()
+ 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
"""
- (J1,J2) = self.factors()
- m = J1.dimension()
- n = J2.dimension()
- V_basis = self.vector_space().basis()
- # Need to specify the dimensions explicitly so that we don't
- # wind up with a zero-by-zero matrix when we want e.g. a
- # two-by-zero matrix (important for composing things).
- I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
- I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
- iota_left = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,self,I1)
- iota_right = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,self,I2)
- return (iota_left, iota_right)
+ 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())
- def inner_product(self, x, y):
- r"""
- The standard Cartesian inner-product.
- We project ``x`` and ``y`` onto our factors, and add up the
- inner-products from the subalgebras.
+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 (FiniteDimensionalEJA,
+ ....: HadamardEJA,
+ ....: JordanSpinEJA,
+ ....: RealSymmetricEJA)
+
+ EXAMPLES:
+
+ This gives us fast characteristic polynomial computations in
+ product algebras, too::
+
+
+ sage: J1 = JordanSpinEJA(2)
+ sage: J2 = RealSymmetricEJA(3)
+ sage: J = cartesian_product([J1,J2])
+ sage: J.characteristic_polynomial_of().degree()
+ 5
+ sage: J.rank()
+ 5
+
+ TESTS:
+
+ The ``cartesian_product()`` function only uses the first factor to
+ decide where the result will live; thus we have to be careful to
+ check that all factors do indeed have a ``rational_algebra()`` method
+ before we construct an algebra that claims to have a rational basis::
+
+ sage: J1 = HadamardEJA(2)
+ sage: jp = lambda X,Y: X*Y
+ sage: ip = lambda X,Y: X[0,0]*Y[0,0]
+ sage: b1 = matrix(QQ, [[1]])
+ sage: J2 = FiniteDimensionalEJA((b1,), jp, ip)
+ sage: cartesian_product([J2,J1]) # factor one not RationalBasisEJA
+ Euclidean Jordan algebra of dimension 1 over Algebraic Real
+ Field (+) Euclidean Jordan algebra of dimension 2 over Algebraic
+ Real Field
+ sage: cartesian_product([J1,J2]) # factor one is RationalBasisEJA
+ Traceback (most recent call last):
+ ...
+ ValueError: factor not a RationalBasisEJA
+
+ """
+ def __init__(self, algebras, **kwargs):
+ if not all( hasattr(r, "rational_algebra") for r in algebras ):
+ raise ValueError("factor not a RationalBasisEJA")
+
+ CartesianProductEJA.__init__(self, algebras, **kwargs)
+
+ @cached_method
+ def rational_algebra(self):
+ if self.base_ring() is QQ:
+ return self
+
+ return cartesian_product([
+ r.rational_algebra() for r in self.cartesian_factors()
+ ])
+
+
+RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
+
+def random_eja(max_dimension=None, *args, **kwargs):
+ r"""
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import random_eja
+
+ TESTS::
+
+ sage: n = ZZ.random_element(1,5)
+ sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False)
+ sage: J.dimension() <= n
+ True
+
+ """
+ # Use the ConcreteEJA default as the total upper bound (regardless
+ # of any whether or not any individual factors set a lower limit).
+ if max_dimension is None:
+ max_dimension = ConcreteEJA._max_random_instance_dimension()
+ J1 = ConcreteEJA.random_instance(max_dimension, *args, **kwargs)
+
+
+ # Roll the dice to see if we attempt a Cartesian product.
+ dice_roll = ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1)
+ new_max_dimension = max_dimension - J1.dimension()
+ if new_max_dimension == 0 or dice_roll != 0:
+ # If it's already as big as we're willing to tolerate, just
+ # return it and don't worry about Cartesian products.
+ return J1
+ else:
+ # Use random_eja() again so we can get more than two factors
+ # if the sub-call also Decides on a cartesian product.
+ J2 = random_eja(new_max_dimension, *args, **kwargs)
+ return cartesian_product([J1,J2])
+
+
+class ComplexSkewSymmetricEJA(RationalBasisEJA, ConcreteEJA):
+ r"""
+ The skew-symmetric EJA of order `2n` described in Faraut and
+ Koranyi's Exercise III.1.b. It has dimension `2n^2 - n`.
+
+ It is (not obviously) isomorphic to the QuaternionHermitianEJA of
+ order `n`, as can be inferred by comparing rank/dimension or
+ explicitly from their "characteristic polynomial of" functions,
+ which just so happen to align nicely.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (ComplexSkewSymmetricEJA,
+ ....: QuaternionHermitianEJA)
+ sage: from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
+
+ EXAMPLES:
+
+ This EJA is isomorphic to the quaternions::
+
+ sage: J = ComplexSkewSymmetricEJA(2, field=QQ, orthonormalize=False)
+ sage: K = QuaternionHermitianEJA(2, field=QQ, orthonormalize=False)
+ sage: jordan_isom_matrix = matrix.diagonal(QQ,[-1,1,1,1,1,-1])
+ sage: phi = FiniteDimensionalEJAOperator(J,K,jordan_isom_matrix)
+ sage: all( phi(x*y) == phi(x)*phi(y)
+ ....: for x in J.gens()
+ ....: for y in J.gens() )
+ True
+ sage: x,y = J.random_elements(2)
+ sage: phi(x*y) == phi(x)*phi(y)
+ True
+
+ TESTS:
+
+ Random elements should satisfy the same conditions that the basis
+ elements do::
+
+ sage: K = ComplexSkewSymmetricEJA.random_instance(field=QQ,
+ ....: orthonormalize=False)
+ sage: x,y = K.random_elements(2)
+ sage: z = x*y
+ sage: x = x.to_matrix()
+ sage: y = y.to_matrix()
+ sage: z = z.to_matrix()
+ sage: all( e.is_skew_symmetric() for e in (x,y,z) )
+ True
+ sage: J = -K.one().to_matrix()
+ sage: all( e*J == J*e.conjugate() for e in (x,y,z) )
+ True
+
+ The power law in Faraut & Koranyi's II.7.a is satisfied.
+ We're in a subalgebra of theirs, but powers are still
+ defined the same::
+
+ sage: K = ComplexSkewSymmetricEJA.random_instance(field=QQ,
+ ....: orthonormalize=False)
+ sage: x = K.random_element()
+ sage: k = ZZ.random_element(5)
+ sage: actual = x^k
+ sage: J = -K.one().to_matrix()
+ sage: expected = K(-J*(J*x.to_matrix())^k)
+ sage: actual == expected
+ True
+
+ """
+ @staticmethod
+ def _max_random_instance_size(max_dimension):
+ # Obtained by solving d = 2n^2 - n, which comes from noticing
+ # that, in 2x2 block form, any element of this algebra has a
+ # free skew-symmetric top-left block, a Hermitian top-right
+ # block, and two bottom blocks that are determined by the top.
+ # The ZZ-int-ZZ thing is just "floor."
+ return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4))
+
+ @classmethod
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
+ """
+ Return a random instance of this type of algebra.
+ """
+ class_max_d = cls._max_random_instance_dimension()
+ if (max_dimension is None or max_dimension > class_max_d):
+ max_dimension = class_max_d
+ max_size = cls._max_random_instance_size(max_dimension)
+ n = ZZ.random_element(max_size + 1)
+ return cls(n, **kwargs)
+
+ @staticmethod
+ def _denormalized_basis(A):
+ """
SETUP::
+ sage: from mjo.hurwitz import ComplexMatrixAlgebra
+ sage: from mjo.eja.eja_algebra import ComplexSkewSymmetricEJA
- sage: from mjo.eja.eja_algebra import (HadamardEJA,
- ....: QuaternionHermitianEJA,
- ....: DirectSumEJA)
-
- EXAMPLE::
-
- sage: J1 = HadamardEJA(3,field=QQ)
- sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
- sage: J = DirectSumEJA(J1,J2)
- sage: x1 = J1.one()
- sage: x2 = x1
- sage: y1 = J2.one()
- sage: y2 = y1
- sage: x1.inner_product(x2)
- 3
- sage: y1.inner_product(y2)
- 2
- sage: J.one().inner_product(J.one())
- 5
+ TESTS:
+
+ The basis elements are all skew-Hermitian::
+
+ sage: d_max = ComplexSkewSymmetricEJA._max_random_instance_dimension()
+ sage: n_max = ComplexSkewSymmetricEJA._max_random_instance_size(d_max)
+ sage: n = ZZ.random_element(n_max + 1)
+ sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ)
+ sage: B = ComplexSkewSymmetricEJA._denormalized_basis(A)
+ sage: all( M.is_skew_symmetric() for M in B)
+ True
+
+ The basis elements ``b`` all satisfy ``b*J == J*b.conjugate()``,
+ as in the definition of the algebra::
+
+ sage: d_max = ComplexSkewSymmetricEJA._max_random_instance_dimension()
+ sage: n_max = ComplexSkewSymmetricEJA._max_random_instance_size(d_max)
+ sage: n = ZZ.random_element(n_max + 1)
+ sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ)
+ sage: I_n = matrix.identity(ZZ, n)
+ sage: J = matrix.block(ZZ, 2, 2, (0, I_n, -I_n, 0), subdivide=False)
+ sage: J = A.from_list(J.rows())
+ sage: B = ComplexSkewSymmetricEJA._denormalized_basis(A)
+ sage: all( b*J == J*b.conjugate() for b in B )
+ True
"""
- (pi_left, pi_right) = self.projections()
- x1 = pi_left(x)
- x2 = pi_right(x)
- y1 = pi_left(y)
- y2 = pi_right(y)
+ es = A.entry_algebra_gens()
+ gen = lambda A,m: A.monomial(m)
+
+ basis = []
+
+ # The size of the blocks. We're going to treat these thing as
+ # 2x2 block matrices,
+ #
+ # [ x1 x2 ]
+ # [ -x2-conj x1-conj ]
+ #
+ # where x1 is skew-symmetric and x2 is Hermitian.
+ #
+ m = A.nrows()/2
+
+ # We only loop through the top half of the matrix, because the
+ # bottom can be constructed from the top.
+ for i in range(m):
+ # First do the top-left block, which is skew-symmetric.
+ # We can compute the bottom-right block in the process.
+ for j in range(i+1):
+ if i != j:
+ # Skew-symmetry implies zeros for (i == j).
+ for e in es:
+ # Top-left block's entry.
+ E_ij = gen(A, (i,j,e))
+ E_ij -= gen(A, (j,i,e))
+
+ # Bottom-right block's entry.
+ F_ij = gen(A, (i+m,j+m,e)).conjugate()
+ F_ij -= gen(A, (j+m,i+m,e)).conjugate()
+
+ basis.append(E_ij + F_ij)
+
+ # Now do the top-right block, which is Hermitian, and compute
+ # the bottom-left block along the way.
+ for j in range(m,i+m+1):
+ if (i+m) == j:
+ # Hermitian matrices have real diagonal entries.
+ # Top-right block's entry.
+ E_ii = gen(A, (i,j,es[0]))
+
+ # Bottom-left block's entry. Don't conjugate
+ # 'cause it's real.
+ E_ii -= gen(A, (i+m,j-m,es[0]))
+ basis.append(E_ii)
+ else:
+ for e in es:
+ # Top-right block's entry. BEWARE! We're not
+ # reflecting across the main diagonal as in
+ # (i,j)~(j,i). We're only reflecting across
+ # the diagonal for the top-right block.
+ E_ij = gen(A, (i,j,e))
+
+ # Shift it back to non-offset coords, transpose,
+ # conjugate, and put it back:
+ #
+ # (i,j) -> (i,j-m) -> (j-m, i) -> (j-m, i+m)
+ E_ij += gen(A, (j-m,i+m,e)).conjugate()
- return (x1.inner_product(y1) + x2.inner_product(y2))
+ # Bottom-left's block's below-diagonal entry.
+ # Just shift the top-right coords down m and
+ # left m.
+ F_ij = -gen(A, (i+m,j-m,e)).conjugate()
+ F_ij += -gen(A, (j,i,e)) # double-conjugate cancels
+ basis.append(E_ij + F_ij)
+ return tuple( basis )
-random_eja = ConcreteEuclideanJordanAlgebra.random_instance
+ @staticmethod
+ @cached_method
+ def _J_matrix(matrix_space):
+ n = matrix_space.nrows() // 2
+ F = matrix_space.base_ring()
+ I_n = matrix.identity(F, n)
+ J = matrix.block(F, 2, 2, (0, I_n, -I_n, 0), subdivide=False)
+ return matrix_space.from_list(J.rows())
+
+ def J_matrix(self):
+ return ComplexSkewSymmetricEJA._J_matrix(self.matrix_space())
+
+ def __init__(self, n, field=AA, **kwargs):
+ # New code; always check the axioms.
+ #if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+ from mjo.hurwitz import ComplexMatrixAlgebra
+ A = ComplexMatrixAlgebra(2*n, scalars=field)
+ J = ComplexSkewSymmetricEJA._J_matrix(A)
+
+ def jordan_product(X,Y):
+ return (X*J*Y + Y*J*X)/2
+
+ def inner_product(X,Y):
+ return (X.conjugate_transpose()*Y).trace().real()
+
+ super().__init__(self._denormalized_basis(A),
+ jordan_product,
+ inner_product,
+ field=field,
+ matrix_space=A,
+ **kwargs)
+
+ # This algebra is conjectured (by me) to be isomorphic to
+ # the quaternion Hermitian EJA of size n, and the rank
+ # would follow from that.
+ #self.rank.set_cache(n)
+ self.one.set_cache( self(-J) )