X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=615b0d494af83a132e1d6343fe05917e4ed2eaf7;hb=2eead94218aad1f63faec9cdeacc30171a880438;hp=a2adbb2bed47499aea53b38fe175547b887fb2a9;hpb=e6bb7aba7677a04c3aaa4d1fdd9dd45185db2067;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index a2adbb2..615b0d4 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 @@ -227,9 +319,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. @@ -1661,13 +1755,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() @@ -1715,25 +1802,35 @@ class ConcreteEJA(FiniteDimensionalEJA): """ @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 ``max_dimension``. If the dimension bound + is omitted, then the ``_max_random_instance_dimension()`` is used + to get a suitable bound. This method should be implemented in each subclass. """ @@ -1937,8 +2034,8 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): 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 @@ -1969,15 +2066,20 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): """ @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, **kwargs): """ Return a random instance of this type of algebra. """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) + 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) return cls(n, **kwargs) def __init__(self, n, field=AA, **kwargs): @@ -2039,8 +2141,8 @@ 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 @@ -2068,6 +2170,7 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): 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 @@ -2087,15 +2190,19 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): 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, **kwargs): """ Return a random instance of this type of algebra. """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) + if max_dimension is None: + max_dimension = cls._max_random_instance_dimension() + n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1) return cls(n, **kwargs) @@ -2125,8 +2232,8 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): 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 @@ -2176,18 +2283,22 @@ class QuaternionHermitianEJA(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 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, **kwargs): """ Return a random instance of this type of algebra. """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) + if max_dimension is None: + max_dimension = cls._max_random_instance_dimension() + n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1) return cls(n, **kwargs) class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): @@ -2278,18 +2389,29 @@ 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, **kwargs): """ Return a random instance of this type of algebra. """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) + if max_dimension is None: + max_dimension = cls._max_random_instance_dimension() + n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1) return cls(n, **kwargs) def __init__(self, n, field=AA, **kwargs): @@ -2410,18 +2532,28 @@ class HadamardEJA(RationalBasisEJA, ConcreteEJA): 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, **kwargs): """ Return a random instance of this type of algebra. """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) + if max_dimension is None: + max_dimension = cls._max_random_instance_dimension() + n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1) return cls(n, **kwargs) @@ -2561,18 +2693,28 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA): 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, **kwargs): """ Return a random instance of this algebra. """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) + if max_dimension is None: + max_dimension = cls._max_random_instance_dimension() + n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1) if n.is_zero(): B = matrix.identity(ZZ, n) return cls(B, **kwargs) @@ -2655,21 +2797,16 @@ class JordanSpinEJA(BilinearFormEJA): # 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, **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) + if max_dimension is None: + max_dimension = cls._max_random_instance_dimension() + n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1) return cls(n, **kwargs) @@ -2944,6 +3081,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)) @@ -3273,15 +3417,23 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA, RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA -def random_eja(*args, **kwargs): - J1 = ConcreteEJA.random_instance(*args, **kwargs) +def random_eja(max_dimension=None, *args, **kwargs): + # 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) - # 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: + + # 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])