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: from mjo.eja.eja_algebra import random_eja
from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
from mjo.eja.eja_utils import _all2list, _mat2vec
+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 FiniteDimensionalEJA(CombinatorialFreeModule):
r"""
A finite-dimensional Euclidean Jordan algebra.
"""
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,
matrix_space=None,
orthonormalize=True,
associative=None,
- cartesian_product=False,
check_field=True,
check_axioms=True,
prefix="b"):
n = len(basis)
if check_field:
- if not field.is_subring(RR):
- # Note: this does return true for the real algebraic
- # field, the rationals, and any quadratic field where
- # we've specified a real embedding.
- raise ValueError("scalar field is not real")
+ 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.
- if not all( jordan_product(bi,bj) == jordan_product(bj,bi)
- for bi in basis
- for bj in basis ):
- raise ValueError("Jordan product is not commutative")
-
- if not all( inner_product(bi,bj) == inner_product(bj,bi)
- for bi in basis
- for bj in basis ):
- raise ValueError("inner-product is not commutative")
-
-
- category = MagmaticAlgebras(field).FiniteDimensional()
- category = category.WithBasis().Unital().Commutative()
+ self._check_input_axioms(basis, jordan_product, inner_product)
if n <= 1:
# All zero- and one-dimensional algebras are just the real
for bj in basis
for bk in basis)
+ category = EuclideanJordanAlgebras(field)
+
if associative:
# Element subalgebras can take advantage of this.
category = category.Associative()
- if cartesian_product:
- # Use join() here because otherwise we only get the
- # "Cartesian product of..." and not the things themselves.
- category = category.join([category,
- category.CartesianProducts()])
# Call the superclass constructor so that we can use its from_vector()
# method to build our multiplication table.
# as well as a subspace W of V spanned by those (vectorized)
# basis elements. The W-coordinates are the coefficients that
# we see in things like x = 1*b1 + 2*b2.
- vector_basis = basis
degree = 0
if n > 0:
# written out as "long vectors."
V = VectorSpace(field, degree)
- # The matrix that will hole the orthonormal -> unorthonormal
- # coordinate transformation.
- self._deortho_matrix = None
+ # 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.
# its own set of non-ambient coordinates (in terms of the
# supplied basis).
vector_basis = tuple( V(_all2list(b)) for b in basis )
- W = V.span_of_basis( vector_basis, check=check_axioms)
+
+ # 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 "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.
+ # Now "self._matrix_span" is the vector space of our
+ # algebra coordinates. The variables "X1", "X2",... 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( U.coordinate_vector(q)
- for q in vector_basis )
+ 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)
- self._multiplication_table = [ [0 for j in range(i+1)]
+ 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
# The jordan product returns a matrixy answer, so we
# have to convert it to the algebra coordinates.
elt = jordan_product(q_i, q_j)
- elt = W.coordinate_vector(V(_all2list(elt)))
+ elt = self._matrix_span.coordinate_vector(V(_all2list(elt)))
self._multiplication_table[i][j] = self.from_vector(elt)
if not orthonormalize:
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.
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 back and forth
Traceback (most recent call last):
...
ValueError: not an element of this algebra
+
"""
msg = "not an element of this algebra"
if elt in self.base_ring():
# that the integer 3 belongs to the space of 2-by-2 matrices.
raise ValueError(msg)
- try:
- # Try to convert a vector into a column-matrix...
+ 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()
- except (AttributeError, TypeError):
- # and ignore failure, because we weren't really expecting
- # a vector as an argument anyway.
- pass
if elt not in self.matrix_space():
raise ValueError(msg)
# is that we're already converting everything to long vectors,
# and that strategy works for tuples as well.
#
- # We pass check=False because the matrix basis is "guaranteed"
- # to be linearly independent... right? Ha ha.
- elt = _all2list(elt)
- V = VectorSpace(self.base_ring(), len(elt))
- W = V.span_of_basis( (V(_all2list(s)) for s in self.matrix_basis()),
- check=False)
+ elt = self._matrix_span.ambient_vector_space()(_all2list(elt))
try:
- coords = W.coordinate_vector(V(elt))
+ coords = self._matrix_span.coordinate_vector(elt)
except ArithmeticError: # vector is not in free module
raise ValueError(msg)
# corresponding to trivial spaces (e.g. it returns only the
# eigenspace corresponding to lambda=1 if you take the
# decomposition relative to the identity element).
- trivial = self.subalgebra(())
+ trivial = self.subalgebra((), check_axioms=False)
J0 = trivial # eigenvalue zero
J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
J1 = trivial # eigenvalue one
a = ( a_i.change_ring(self.base_ring())
for a_i in self._rational_algebra._charpoly_coefficients() )
- if self._deortho_matrix is None:
- # This can happen if our base ring was, say, AA and we
- # chose not to (or didn't need to) orthonormalize. It's
- # still faster to do the computations over QQ even if
- # the numbers in the boxes stay the same.
- return tuple(a)
-
# Otherwise, convert the coordinate variables back to the
# deorthonormalized ones.
R = self.coordinate_polynomial_ring()
def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
Return a random instance of this type of algebra whose dimension
- is less than or equal to ``max_dimension``. If the dimension bound
- is omitted, then the ``_max_random_instance_dimension()`` is used
- to get a suitable bound.
+ 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
+
"""
from sage.misc.prandom import choice
eja_class = choice(cls.__subclasses__())
# These all bubble up to the RationalBasisEJA superclass
# constructor, so any (kw)args valid there are also valid
# here.
- return eja_class.random_instance(*args, **kwargs)
+ return eja_class.random_instance(max_dimension, *args, **kwargs)
class MatrixEJA(FiniteDimensionalEJA):
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
@staticmethod
def _max_random_instance_size(max_dimension):
# Obtained by solving d = (n^2 + n)/2.
- return int(ZZ(8*max_dimension + 1).sqrt()/2 - 1/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, max_dimension=None, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
Return a random instance of this type of algebra.
"""
- if max_dimension is None:
- max_dimension = cls._max_random_instance_dimension()
- max_size = cls._max_random_instance_size(max_dimension) + 1
- n = ZZ.random_element(max_size)
+ 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):
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
@staticmethod
def _max_random_instance_size(max_dimension):
# Obtained by solving d = n^2.
- return int(ZZ(max_dimension).sqrt())
+ # The ZZ-int-ZZ thing is just "floor."
+ return ZZ(int(ZZ(max_dimension).sqrt()))
@classmethod
- def random_instance(cls, max_dimension=None, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
Return a random instance of this type of algebra.
"""
- if max_dimension is None:
- max_dimension = cls._max_random_instance_dimension()
- n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 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)
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 maximum rank of a random QuaternionHermitianEJA.
"""
# Obtained by solving d = 2n^2 - n.
- return int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4)
+ # 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, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
Return a random instance of this type of algebra.
"""
- if max_dimension is None:
- max_dimension = cls._max_random_instance_dimension()
- n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 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 OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
return 0
@classmethod
- def random_instance(cls, max_dimension=None, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
Return a random instance of this type of algebra.
"""
- if max_dimension is None:
- max_dimension = cls._max_random_instance_dimension()
- n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 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):
return max_dimension
@classmethod
- def random_instance(cls, max_dimension=None, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
Return a random instance of this type of algebra.
"""
- if max_dimension is None:
- max_dimension = cls._max_random_instance_dimension()
- n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 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)
return max_dimension
@classmethod
- def random_instance(cls, max_dimension=None, **kwargs):
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
"""
Return a random instance of this algebra.
"""
- if max_dimension is None:
- max_dimension = cls._max_random_instance_dimension()
- n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 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
super().__init__(B, *args, **kwargs)
@classmethod
- def random_instance(cls, max_dimension=None, **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``.
"""
- if max_dimension is None:
- max_dimension = cls._max_random_instance_dimension()
- n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 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)
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)
sage: from mjo.eja.eja_algebra import (random_eja,
....: CartesianProductEJA,
+ ....: ComplexHermitianEJA,
....: HadamardEJA,
....: JordanSpinEJA,
....: RealSymmetricEJA)
| 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:
All factors must share the same base field::
sage: expected = J.one() # long time
sage: actual == expected # long time
True
-
"""
- Element = FiniteDimensionalEJAElement
-
-
def __init__(self, factors, **kwargs):
m = len(factors)
if m == 0:
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 )
-
- # Compute my matrix space. This category isn't perfect, but
- # is good enough for what we need to do.
+ 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
- MS = CartesianProduct(MS_factors, MS_cat)
+ self._matrix_space = CartesianProduct(MS_factors, MS_cat)
- basis = []
- zero = MS.zero()
+ 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
- basis.append(z)
+ self._matrix_basis.append(z)
- basis = tuple( MS(b) for b in basis )
+ self._matrix_basis = tuple( self._matrix_space(b)
+ for b in self._matrix_basis )
+ n = len(self._matrix_basis)
- # Define jordan/inner products that operate on that matrix_basis.
- def jordan_product(x,y):
- return MS(tuple(
- (factors[i](x[i])*factors[i](y[i])).to_matrix()
- for i in range(m)
- ))
-
- def inner_product(x, y):
- return sum(
- factors[i](x[i]).inner_product(factors[i](y[i]))
- for i in range(m)
- )
+ # We already have what we need for the super-superclass constructor.
+ CombinatorialFreeModule.__init__(self,
+ field,
+ range(n),
+ prefix="b",
+ category=category,
+ bracket=False)
- # There's no need to check the field since it already came
- # from an EJA. Likewise the axioms are guaranteed to be
- # satisfied, unless the guy writing this class sucks.
- #
- # If you want the basis to be orthonormalized, orthonormalize
- # the factors.
- FiniteDimensionalEJA.__init__(self,
- basis,
- jordan_product,
- inner_product,
- field=field,
- matrix_space=MS,
- orthonormalize=False,
- associative=associative,
- cartesian_product=True,
- check_field=False,
- check_axioms=False)
+ # 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]
+ )
+
+ self._inner_product_matrix = matrix.block_diagonal(
+ [J._inner_product_matrix for J in factors]
+ )
+
+ # 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)
return cartesian_product.symbol.join("%s" % factor
for factor in self._sets)
- def matrix_space(self):
- r"""
- Return the space that our matrix basis lives in as a Cartesian
- product.
-
- We don't simply use the ``cartesian_product()`` functor here
- because it acts differently on SageMath MatrixSpaces and our
- custom MatrixAlgebras, which are CombinatorialFreeModules. We
- always want the result to be represented (and indexed) as
- an ordered tuple.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
- ....: HadamardEJA,
- ....: OctonionHermitianEJA,
- ....: RealSymmetricEJA)
-
- EXAMPLES::
-
- sage: J1 = HadamardEJA(1)
- sage: J2 = RealSymmetricEJA(2)
- sage: J = cartesian_product([J1,J2])
- sage: J.matrix_space()
- The Cartesian product of (Full MatrixSpace of 1 by 1 dense
- matrices over Algebraic Real Field, Full MatrixSpace of 2
- by 2 dense matrices over Algebraic Real Field)
-
- ::
-
- sage: J1 = ComplexHermitianEJA(1)
- sage: J2 = ComplexHermitianEJA(1)
- sage: J = cartesian_product([J1,J2])
- sage: J.one().to_matrix()[0]
- +---+
- | 1 |
- +---+
- sage: J.one().to_matrix()[1]
- +---+
- | 1 |
- +---+
-
- ::
-
- sage: J1 = OctonionHermitianEJA(1)
- sage: J2 = OctonionHermitianEJA(1)
- sage: J = cartesian_product([J1,J2])
- sage: J.one().to_matrix()[0]
- +----+
- | e0 |
- +----+
- sage: J.one().to_matrix()[1]
- +----+
- | e0 |
- +----+
-
- """
- return super().matrix_space()
-
@cached_method
def cartesian_projection(self, i):
RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
def random_eja(max_dimension=None, *args, **kwargs):
+ r"""
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import random_eja
+
+ TESTS::
+
+ sage: set_random_seed()
+ 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: