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
# 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
# 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)
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()
"""
@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
+
"""
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():
- 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, field=AA, **kwargs):
A = MatrixSpace(field, n)
super().__init__(A, **kwargs)
+ from mjo.eja.eja_cache import real_symmetric_eja_coeffs
+ a = real_symmetric_eja_coeffs(self)
+ if a is not None:
+ if self._rational_algebra is None:
+ self._charpoly_coefficients.set_cache(a)
+ else:
+ self._rational_algebra._charpoly_coefficients.set_cache(a)
+
class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
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
sage: ComplexHermitianEJA(0)
Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+
"""
def __init__(self, n, field=AA, **kwargs):
# We know this is a valid EJA, but will double-check
A = ComplexMatrixAlgebra(n, scalars=field)
super().__init__(A, **kwargs)
+ from mjo.eja.eja_cache import complex_hermitian_eja_coeffs
+ a = complex_hermitian_eja_coeffs(self)
+ if a is not None:
+ if self._rational_algebra is None:
+ self._charpoly_coefficients.set_cache(a)
+ else:
+ 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)
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
@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
+ # 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))
@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 OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
"""
@staticmethod
- def _max_random_instance_size():
+ def _max_random_instance_size(max_dimension):
r"""
The maximum rank of a random QuaternionHermitianEJA.
"""
- return 1 # Dimension 1
+ # 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):
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)
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
# can pass in a field!
super().__init__(B, *args, **kwargs)
- @staticmethod
- def _max_random_instance_size():
- r"""
- The maximum dimension of a random JordanSpinEJA.
- """
- return 5
-
@classmethod
- def random_instance(cls, **kwargs):
+ 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)
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)
check_field=False,
check_axioms=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.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))
RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
-def random_eja(*args, **kwargs):
- J1 = ConcreteEJA.random_instance(*args, **kwargs)
+def random_eja(max_dimension=None, *args, **kwargs):
+ r"""
- # This might make Cartesian products appear roughly as often as
- # any other ConcreteEJA.
- if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0:
- # Use random_eja() again so we can get more than two factors.
- J2 = random_eja(*args, **kwargs)
- J = cartesian_product([J1,J2])
- return J
- else:
+ 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:
+ 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])
projects / sage.d.git / blobdiff