X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=4a1c6f94f589c80b9021b71ee01693765aabcaeb;hb=82d3e5bb067d73eb0aa557d234da0f9723b456e8;hp=db1f4769ad84068903356bc119749ad90a64dc88;hpb=046d2b5d664f7b6794e81e7ebf0fb224c0c3d52c;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index db1f476..4a1c6f9 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -52,6 +52,98 @@ 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: from mjo.eja.eja_algebra import random_eja @@ -217,7 +309,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # 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: @@ -227,9 +318,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # 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. @@ -254,17 +347,22 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # 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 @@ -285,7 +383,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # 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: @@ -592,8 +690,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): 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. @@ -632,6 +730,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): 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 @@ -656,6 +764,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): 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(): @@ -665,13 +774,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # 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, 'column'): + # Convert a vector into a column-matrix... elt = elt.column() - except (AttributeError, TypeError): - # and ignore failure, because we weren't really expecting - # a vector as an argument anyway. - pass if elt not in self.matrix_space(): raise ValueError(msg) @@ -688,15 +798,10 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # 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) @@ -1661,13 +1766,6 @@ class RationalBasisEJA(FiniteDimensionalEJA): 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() @@ -1741,11 +1839,36 @@ class ConcreteEJA(FiniteDimensionalEJA): 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__()) @@ -1753,7 +1876,7 @@ class ConcreteEJA(FiniteDimensionalEJA): # 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): @@ -1985,14 +2108,15 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): 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): @@ -2109,13 +2233,15 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): 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) @@ -2205,13 +2331,15 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): 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): @@ -2318,13 +2446,15 @@ 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): @@ -2460,13 +2590,15 @@ class HadamardEJA(RationalBasisEJA, ConcreteEJA): 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) @@ -2621,13 +2753,16 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA): 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) @@ -2638,6 +2773,7 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA): 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 @@ -2711,15 +2847,17 @@ class JordanSpinEJA(BilinearFormEJA): 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) @@ -2776,9 +2914,12 @@ class TrivialEJA(RationalBasisEJA, ConcreteEJA): 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=`` argument, we don't try to pass + # it on to the algebra constructor. return cls(**kwargs) @@ -2994,6 +3135,13 @@ class CartesianProductEJA(FiniteDimensionalEJA): 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)) @@ -3324,6 +3472,21 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA, 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: