X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=8ca501d6ba944a539901d384a78a5c8a0f45a7ca;hb=3f8818f3eca0b2ea388130ff805875012cf902cb;hp=4091d03fd6a40351e7d4c2789847d1db6fa0e807;hpb=4d0f89a814fe6ab91a17af023c35caefaada2893;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 4091d03..8ca501d 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -57,11 +57,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): prefix='e', category=None, natural_basis=None, - check=True): + check_field=True, + check_axioms=True): """ SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja) + sage: from mjo.eja.eja_algebra import ( + ....: FiniteDimensionalEuclideanJordanAlgebra, + ....: JordanSpinEJA, + ....: random_eja) EXAMPLES: @@ -75,20 +79,33 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): TESTS: - The ``field`` we're given must be real:: + The ``field`` we're given must be real with ``check_field=True``:: sage: JordanSpinEJA(2,QQbar) Traceback (most recent call last): ... - ValueError: field is not real + ValueError: scalar field is not real + + The multiplication table must be square with ``check_axioms=True``:: + + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),())) + Traceback (most recent call last): + ... + ValueError: multiplication table is not square """ - if check: + if check_field: if not field.is_subring(RR): # Note: this does return true for the real algebraic - # field, and any quadratic field where we've specified - # a real embedding. - raise ValueError('field is not real') + # field, the rationals, and any quadratic field where + # we've specified a real embedding. + raise ValueError("scalar field is not real") + + # The multiplication table had better be square + n = len(mult_table) + if check_axioms: + if not all( len(l) == n for l in mult_table ): + raise ValueError("multiplication table is not square") self._natural_basis = natural_basis @@ -98,7 +115,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): fda = super(FiniteDimensionalEuclideanJordanAlgebra, self) fda.__init__(field, - range(len(mult_table)), + range(n), prefix=prefix, category=category) self.print_options(bracket='') @@ -114,6 +131,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): for ls in mult_table ] + if check_axioms: + if not self._is_commutative(): + raise ValueError("algebra is not commutative") + 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 _element_constructor_(self, elt): """ @@ -193,7 +217,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return self.from_vector(coords) @staticmethod - def _max_test_case_size(): + def _max_random_instance_size(): """ Return an integer "size" that is an upper bound on the size of this algebra when it is used in a random test @@ -209,7 +233,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): interpreted to be far less than the dimension) should override with a smaller number. """ - return 5 + raise NotImplementedError def _repr_(self): """ @@ -235,11 +259,74 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): def product_on_basis(self, i, j): return self._multiplication_table[i][j] + 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, 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(self): + def characteristic_polynomial_of(self): """ - Return a characteristic polynomial that works for all elements - of this algebra. + Return the algebra's "characteristic polynomial of" function, + which is itself a multivariate polynomial that, when evaluated + at the coordinates of some algebra element, returns that + element's characteristic polynomial. The resulting polynomial has `n+1` variables, where `n` is the dimension of this algebra. The first `n` variables correspond to @@ -259,7 +346,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Alizadeh, Example 11.11:: sage: J = JordanSpinEJA(3) - sage: p = J.characteristic_polynomial(); p + sage: p = J.characteristic_polynomial_of(); p X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2 sage: xvec = J.one().to_vector() sage: p(*xvec) @@ -272,7 +359,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): any argument:: sage: J = TrivialEJA() - sage: J.characteristic_polynomial() + sage: J.characteristic_polynomial_of() 1 """ @@ -438,8 +525,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): """ Return the matrix space in which this algebra's natural basis elements live. + + Generally this will be an `n`-by-`1` column-vector space, + except when the algebra is trivial. There it's `n`-by-`n` + (where `n` is zero), to ensure that two elements of the + natural basis space (empty matrices) can be multiplied. """ - if self._natural_basis is None or len(self._natural_basis) == 0: + if self.is_trivial(): + return MatrixSpace(self.base_ring(), 0) + elif self._natural_basis is None or len(self._natural_basis) == 0: return MatrixSpace(self.base_ring(), self.dimension(), 1) else: return self._natural_basis[0].matrix_space() @@ -505,19 +599,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # appeal to the "long vectors" isometry. oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ] - # Now we use basis linear algebra to find the coefficients, + # Now we use basic linear algebra to find the coefficients, # of the matrices-as-vectors-linear-combination, which should # work for the original algebra basis too. - A = matrix.column(self.base_ring(), oper_vecs) + A = matrix(self.base_ring(), oper_vecs) # 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()) ) - # Now if there's an identity element in the algebra, this should work. - coeffs = A.solve_right(b) - return self.linear_combination(zip(self.gens(), coeffs)) + # Now if there's an identity element in the algebra, this + # should work. We solve on the left to avoid having to + # transpose the matrix "A". + return self.from_vector(A.solve_left(b)) def peirce_decomposition(self, c): @@ -572,6 +667,25 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Vector space of degree 6 and dimension 2... sage: J1 Euclidean Jordan algebra of dimension 3... + sage: J0.one().natural_representation() + [0 0 0] + [0 0 0] + [0 0 1] + sage: orig_df = AA.options.display_format + sage: AA.options.display_format = 'radical' + sage: J.from_vector(J5.basis()[0]).natural_representation() + [ 0 0 1/2*sqrt(2)] + [ 0 0 0] + [1/2*sqrt(2) 0 0] + sage: J.from_vector(J5.basis()[1]).natural_representation() + [ 0 0 0] + [ 0 0 1/2*sqrt(2)] + [ 0 1/2*sqrt(2) 0] + sage: AA.options.display_format = orig_df + sage: J1.one().natural_representation() + [1 0 0] + [0 1 0] + [0 0 0] TESTS: @@ -586,9 +700,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J1.superalgebra() == J and J1.dimension() == J.dimension() True - The identity elements in the two subalgebras are the - projections onto their respective subspaces of the - superalgebra's identity element:: + The decomposition is into eigenspaces, and its components are + therefore necessarily orthogonal. Moreover, the identity + 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() @@ -598,6 +713,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): ....: x = J.random_element() sage: c = x.subalgebra_idempotent() sage: J0,J5,J1 = J.peirce_decomposition(c) + sage: ipsum = 0 + sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()): + ....: w = w.superalgebra_element() + ....: y = J.from_vector(y) + ....: z = z.superalgebra_element() + ....: ipsum += w.inner_product(y).abs() + ....: ipsum += w.inner_product(z).abs() + ....: ipsum += y.inner_product(z).abs() + sage: ipsum + 0 sage: J1(c) == J1.one() True sage: J0(J.one() - c) == J0.one() @@ -622,7 +747,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): J5 = eigspace else: gens = tuple( self.from_vector(b) for b in eigspace.basis() ) - subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, gens) + subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, + gens, + check_axioms=False) if eigval == 0: J0 = subalg elif eigval == 1: @@ -633,10 +760,61 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return (J0, J5, J1) - def random_elements(self, count): + def random_element(self, thorough=False): + r""" + Return a random element of this algebra. + + Our algebra superclass method only returns a linear + combination of at most two basis elements. We instead + want the vector space "random element" method that + returns a more diverse selection. + + INPUT: + + - ``thorough`` -- (boolean; default False) whether or not we + should generate irrational coefficients for the random + element when our base ring is irrational; this slows the + algebra operations to a crawl, but any truly random method + should include them + + """ + # For a general base ring... maybe we can trust this to do the + # right thing? Unlikely, but. + V = self.vector_space() + v = V.random_element() + + if self.base_ring() is AA: + # The "random element" method of the algebraic reals is + # stupid at the moment, and only returns integers between + # -2 and 2, inclusive: + # + # https://trac.sagemath.org/ticket/30875 + # + # Instead, we implement our own "random vector" method, + # and then coerce that into the algebra. We use the vector + # space degree here instead of the dimension because a + # subalgebra could (for example) be spanned by only two + # vectors, each with five coordinates. We need to + # generate all five coordinates. + if thorough: + v *= QQbar.random_element().real() + else: + v *= QQ.random_element() + + return self.from_vector(V.coordinate_vector(v)) + + def random_elements(self, count, thorough=False): """ Return ``count`` random elements as a tuple. + INPUT: + + - ``thorough`` -- (boolean; default False) whether or not we + should generate irrational coefficients for the random + elements when our base ring is irrational; this slows the + algebra operations to a crawl, but any truly random method + should include them + SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA @@ -651,7 +829,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): True """ - return tuple( self.random_element() for idx in range(count) ) + return tuple( self.random_element(thorough) + for idx in range(count) ) @classmethod def random_instance(cls, field=AA, **kwargs): @@ -661,12 +840,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Beware, this will crash for "most instances" because the constructor below looks wrong. """ - if cls is TrivialEJA: - # The TrivialEJA class doesn't take an "n" argument because - # there's only one. - return cls(field) - - n = ZZ.random_element(cls._max_test_case_size()) + 1 + n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, field, **kwargs) @cached_method @@ -829,81 +1003,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Element = FiniteDimensionalEuclideanJordanAlgebraElement -class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra): - """ - Return the Euclidean Jordan Algebra corresponding to the set - `R^n` under the Hadamard 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. - - SETUP:: - - sage: from mjo.eja.eja_algebra import HadamardEJA - - EXAMPLES: - - 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 - 0 - sage: e0*e2 - 0 - sage: e1*e1 - e1 - sage: e1*e2 - 0 - sage: e2*e2 - e2 - - TESTS: - - We can change the generator prefix:: - - sage: HadamardEJA(3, prefix='r').gens() - (r0, r1, r2) - - """ - def __init__(self, n, field=AA, **kwargs): - V = VectorSpace(field, n) - mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ] - for i in range(n) ] - - fdeja = super(HadamardEJA, self) - fdeja.__init__(field, mult_table, **kwargs) - self.rank.set_cache(n) - - def inner_product(self, x, y): - """ - Faster to reimplement than to use natural representations. - - SETUP:: - - sage: from mjo.eja.eja_algebra import HadamardEJA - - TESTS: - 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: X = x.natural_representation() - sage: Y = y.natural_representation() - sage: x.inner_product(y) == J.natural_inner_product(X,Y) - True - - """ - return x.to_vector().inner_product(y.to_vector()) - - -def random_eja(field=AA, nontrivial=False): +def random_eja(field=AA): """ Return a "random" finite-dimensional Euclidean Jordan Algebra. @@ -917,24 +1018,79 @@ def random_eja(field=AA, nontrivial=False): Euclidean Jordan algebra of dimension... """ - eja_classes = [HadamardEJA, - JordanSpinEJA, - RealSymmetricEJA, - ComplexHermitianEJA, - QuaternionHermitianEJA] - if not nontrivial: - eja_classes.append(TrivialEJA) - classname = choice(eja_classes) + classname = choice([TrivialEJA, + HadamardEJA, + JordanSpinEJA, + RealSymmetricEJA, + ComplexHermitianEJA, + QuaternionHermitianEJA]) return classname.random_instance(field=field) +class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): + r""" + Algebras whose basis consists of vectors with rational + entries. Equivalently, algebras whose multiplication tables + contain only rational coefficients. + + When an EJA has a basis that can be made rational, we can speed up + the computation of its characteristic polynomial by doing it over + ``QQ``. All of the named EJA constructors that we provide fall + into this category. + """ + @cached_method + def _charpoly_coefficients(self): + r""" + Override the parent method with something that tries to compute + over a faster (non-extension) field. + + SETUP:: + + sage: from mjo.eja.eja_algebra import JordanSpinEJA + + EXAMPLES: + + The base ring of the resulting polynomial coefficients is what + it should be, and not the rationals (unless the algebra was + already over the rationals):: + + sage: J = JordanSpinEJA(3) + sage: J._charpoly_coefficients() + (X1^2 - X2^2 - X3^2, -2*X1) + sage: a0 = J._charpoly_coefficients()[0] + sage: J.base_ring() + Algebraic Real Field + sage: a0.base_ring() + Algebraic Real Field + + """ + if self.base_ring() is QQ: + # There's no need to construct *another* algebra over the + # rationals if this one is already over the rationals. + superclass = super(RationalBasisEuclideanJordanAlgebra, self) + return superclass._charpoly_coefficients() + + mult_table = tuple( + map(lambda x: x.to_vector(), ls) + for ls in self._multiplication_table + ) + + # Do the computation over the rationals. The answer will be + # the same, because our basis coordinates are (essentially) + # rational. + J = FiniteDimensionalEuclideanJordanAlgebra(QQ, + mult_table, + check_field=False, + check_axioms=False) + a = J._charpoly_coefficients() + return tuple(map(lambda x: x.change_ring(self.base_ring()), a)) class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): @staticmethod - def _max_test_case_size(): + def _max_random_instance_size(): # Play it safe, since this will be squared and the underlying # field can have dimension 4 (quaternions) too. return 2 @@ -968,9 +1124,10 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): Qs = self.multiplication_table_from_matrix_basis(basis) - fdeja = super(MatrixEuclideanJordanAlgebra, self) - fdeja.__init__(field, Qs, natural_basis=basis, **kwargs) - return + super(MatrixEuclideanJordanAlgebra, self).__init__(field, + Qs, + natural_basis=basis, + **kwargs) @cached_method @@ -979,38 +1136,44 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): Override the parent method with something that tries to compute over a faster (non-extension) field. """ - if self._basis_normalizers is None: - # We didn't normalize, so assume that the basis we started - # with had entries in a nice field. + if self._basis_normalizers is None or self.base_ring() is QQ: + # We didn't normalize, or the basis we started with had + # entries in a nice field already. Just compute the thing. return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients() - else: - basis = ( (b/n) for (b,n) in zip(self.natural_basis(), - self._basis_normalizers) ) - - # Do this over the rationals and convert back at the end. - # Only works because we know the entries of the basis are - # integers. - J = MatrixEuclideanJordanAlgebra(QQ, - basis, - normalize_basis=False) - a = J._charpoly_coefficients() - - # Unfortunately, changing the basis does change the - # coefficients of the characteristic polynomial, but since - # these are really the coefficients of the "characteristic - # polynomial of" function, everything is still nice and - # unevaluated. It's therefore "obvious" how scaling the - # basis affects the coordinate variables X1, X2, et - # cetera. Scaling the first basis vector up by "n" adds a - # factor of 1/n into every "X1" term, for example. So here - # we simply undo the basis_normalizer scaling that we - # performed earlier. - # - # TODO: make this access safe. - XS = a[0].variables() - subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i] - for i in range(len(XS)) } - return tuple( a_i.subs(subs_dict) for a_i in a ) + + basis = ( (b/n) for (b,n) in zip(self.natural_basis(), + self._basis_normalizers) ) + + # Do this over the rationals and convert back at the end. + # Only works because we know the entries of the basis are + # integers. The argument ``check_axioms=False`` is required + # because the trace inner-product method for this + # class is a stub and can't actually be checked. + J = MatrixEuclideanJordanAlgebra(QQ, + basis, + normalize_basis=False, + check_field=False, + check_axioms=False) + a = J._charpoly_coefficients() + + # Unfortunately, changing the basis does change the + # coefficients of the characteristic polynomial, but since + # these are really the coefficients of the "characteristic + # polynomial of" function, everything is still nice and + # unevaluated. It's therefore "obvious" how scaling the + # basis affects the coordinate variables X1, X2, et + # cetera. Scaling the first basis vector up by "n" adds a + # factor of 1/n into every "X1" term, for example. So here + # we simply undo the basis_normalizer scaling that we + # performed earlier. + # + # The a[0] access here is safe because trivial algebras + # won't have any basis normalizers and therefore won't + # make it to this "else" branch. + XS = a[0].parent().gens() + subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i] + for i in range(len(XS)) } + return tuple( a_i.subs(subs_dict) for a_i in a ) @staticmethod @@ -1028,6 +1191,9 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): # is supposed to hold the entire long vector, and the subspace W # of V will be spanned by the vectors that arise from symmetric # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3. + if len(basis) == 0: + return [] + field = basis[0].base_ring() dimension = basis[0].nrows() @@ -1074,16 +1240,11 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): Yu = cls.real_unembed(Y) tr = (Xu*Yu).trace() - if tr in RLF: - # It's real already. - return tr - - # Otherwise, try the thing that works for complex numbers; and - # if that doesn't work, the thing that works for quaternions. try: - return tr.vector()[0] # real part, imag part is index 1 + # Works in QQ, AA, RDF, et cetera. + return tr.real() except AttributeError: - # A quaternions doesn't have a vector() method, but does + # 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] @@ -1141,7 +1302,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra): The dimension of this algebra is `(n^2 + n) / 2`:: sage: set_random_seed() - sage: n_max = RealSymmetricEJA._max_test_case_size() + sage: n_max = RealSymmetricEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) sage: J = RealSymmetricEJA(n) sage: J.dimension() == (n^2 + n)/2 @@ -1185,6 +1346,11 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra): sage: x.operator().matrix().is_symmetric() True + We can construct the (trivial) algebra of rank zero:: + + sage: RealSymmetricEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ @classmethod def _denormalized_basis(cls, n, field): @@ -1219,13 +1385,16 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra): @staticmethod - def _max_test_case_size(): + def _max_random_instance_size(): return 4 # Dimension 10 def __init__(self, n, field=AA, **kwargs): basis = self._denormalized_basis(n, field) - super(RealSymmetricEJA, self).__init__(field, basis, **kwargs) + super(RealSymmetricEJA, self).__init__(field, + basis, + check_axioms=False, + **kwargs) self.rank.set_cache(n) @@ -1262,7 +1431,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): Embedding is a homomorphism (isomorphism, in fact):: sage: set_random_seed() - sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size() + sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_random_instance_size() sage: n = ZZ.random_element(n_max) sage: F = QuadraticField(-1, 'I') sage: X = random_matrix(F, n) @@ -1414,7 +1583,7 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): The dimension of this algebra is `n^2`:: sage: set_random_seed() - sage: n_max = ComplexHermitianEJA._max_test_case_size() + sage: n_max = ComplexHermitianEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) sage: J = ComplexHermitianEJA(n) sage: J.dimension() == n^2 @@ -1458,6 +1627,11 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): sage: x.operator().matrix().is_symmetric() True + We can construct the (trivial) algebra of rank zero:: + + sage: ComplexHermitianEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ @classmethod @@ -1516,7 +1690,10 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): def __init__(self, n, field=AA, **kwargs): basis = self._denormalized_basis(n,field) - super(ComplexHermitianEJA,self).__init__(field, basis, **kwargs) + super(ComplexHermitianEJA,self).__init__(field, + basis, + check_axioms=False, + **kwargs) self.rank.set_cache(n) @@ -1550,7 +1727,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): Embedding is a homomorphism (isomorphism, in fact):: sage: set_random_seed() - sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size() + sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_random_instance_size() sage: n = ZZ.random_element(n_max) sage: Q = QuaternionAlgebra(QQ,-1,-1) sage: X = random_matrix(Q, n) @@ -1709,7 +1886,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): The dimension of this algebra is `2*n^2 - n`:: sage: set_random_seed() - sage: n_max = QuaternionHermitianEJA._max_test_case_size() + sage: n_max = QuaternionHermitianEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) sage: J = QuaternionHermitianEJA(n) sage: J.dimension() == 2*(n^2) - n @@ -1753,6 +1930,11 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): sage: x.operator().matrix().is_symmetric() True + We can construct the (trivial) algebra of rank zero:: + + sage: QuaternionHermitianEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ @classmethod def _denormalized_basis(cls, n, field): @@ -1812,18 +1994,106 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): def __init__(self, n, field=AA, **kwargs): basis = self._denormalized_basis(n,field) - super(QuaternionHermitianEJA,self).__init__(field, basis, **kwargs) + super(QuaternionHermitianEJA,self).__init__(field, + basis, + check_axioms=False, + **kwargs) self.rank.set_cache(n) -class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra): +class HadamardEJA(RationalBasisEuclideanJordanAlgebra): + """ + Return the Euclidean Jordan Algebra corresponding to the set + `R^n` under the Hadamard 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. + + SETUP:: + + sage: from mjo.eja.eja_algebra import HadamardEJA + + EXAMPLES: + + 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 + 0 + sage: e0*e2 + 0 + sage: e1*e1 + e1 + sage: e1*e2 + 0 + sage: e2*e2 + e2 + + TESTS: + + We can change the generator prefix:: + + sage: HadamardEJA(3, prefix='r').gens() + (r0, r1, r2) + + """ + def __init__(self, n, field=AA, **kwargs): + V = VectorSpace(field, n) + mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ] + for i in range(n) ] + + super(HadamardEJA, self).__init__(field, + mult_table, + check_axioms=False, + **kwargs) + self.rank.set_cache(n) + + @staticmethod + def _max_random_instance_size(): + return 5 + + def inner_product(self, x, y): + """ + Faster to reimplement than to use natural representations. + + SETUP:: + + sage: from mjo.eja.eja_algebra import HadamardEJA + + TESTS: + + 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: X = x.natural_representation() + sage: Y = y.natural_representation() + sage: x.inner_product(y) == J.natural_inner_product(X,Y) + True + + """ + return x.to_vector().inner_product(y.to_vector()) + + +class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): 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 = - (x0*y0 + , x0*y_bar + y0*x_bar)`` where ``B`` is a - symmetric positive-definite "bilinear form" matrix. It has - dimension `n` over the reals, and reduces to the ``JordanSpinEJA`` - when ``B`` is the identity matrix of order ``n-1``. + (,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is + a symmetric positive-definite "bilinear form" matrix. Its + dimension is the size of `B`, and it has rank two in dimensions + larger than two. It reduces to the ``JordanSpinEJA`` when `B` is + the identity matrix of order ``n``. + + We insist that the one-by-one upper-left identity block of `B` be + passed in as well so that we can be passed a matrix of size zero + to construct a trivial algebra. SETUP:: @@ -1835,7 +2105,8 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra): When no bilinear form is specified, the identity matrix is used, and the resulting algebra is the Jordan spin algebra:: - sage: J0 = BilinearFormEJA(3) + sage: B = matrix.identity(AA,3) + sage: J0 = BilinearFormEJA(B) sage: J1 = JordanSpinEJA(3) sage: J0.multiplication_table() == J0.multiplication_table() True @@ -1844,7 +2115,8 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra): We can create a zero-dimensional algebra:: - sage: J = BilinearFormEJA(0) + sage: B = matrix.identity(AA,0) + sage: J = BilinearFormEJA(B) sage: J.basis() Finite family {} @@ -1856,8 +2128,11 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: set_random_seed() sage: n = ZZ.random_element(5) sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular') - sage: B = M.transpose()*M - sage: J = BilinearFormEJA(n, B=B) + sage: B11 = matrix.identity(QQ,1) + sage: B22 = M.transpose()*M + sage: B = block_matrix(2,2,[ [B11,0 ], + ....: [0, B22 ] ]) + sage: J = BilinearFormEJA(B) sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis() sage: V = J.vector_space() sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list())) @@ -1871,11 +2146,12 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: actual == expected True """ - def __init__(self, n, field=AA, B=None, **kwargs): - if B is None: - self._B = matrix.identity(field, max(0,n-1)) - else: - self._B = B + def __init__(self, B, field=AA, **kwargs): + self._B = B + n = B.nrows() + + if not B.is_positive_definite(): + raise TypeError("matrix B is not positive-definite") V = VectorSpace(field, n) mult_table = [[V.zero() for j in range(n)] for i in range(n)] @@ -1887,7 +2163,7 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra): xbar = x[1:] y0 = y[0] ybar = y[1:] - z0 = x0*y0 + (self._B*xbar).inner_product(ybar) + z0 = (B*x).inner_product(y) zbar = y0*xbar + x0*ybar z = V([z0] + zbar.list()) mult_table[i][j] = z @@ -1895,10 +2171,41 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra): # 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). - fdeja = super(BilinearFormEJA, self) - fdeja.__init__(field, mult_table, **kwargs) + super(BilinearFormEJA, self).__init__(field, + mult_table, + check_axioms=False, + **kwargs) self.rank.set_cache(min(n,2)) + @staticmethod + def _max_random_instance_size(): + return 5 + + @classmethod + def random_instance(cls, field=AA, **kwargs): + """ + Return a random instance of this algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + if n == 0: + # Special case needed since we use (n-1) below. + B = matrix.identity(field, 0) + return cls(B, field, **kwargs) + + B11 = matrix.identity(field,1) + M = matrix.random(field, n-1) + I = matrix.identity(field, n-1) + alpha = field.zero() + while alpha.is_zero(): + alpha = field.random_element().abs() + B22 = M.transpose()*M + alpha*I + + from sage.matrix.special import block_matrix + B = block_matrix(2,2, [ [B11, ZZ(0) ], + [ZZ(0), B22 ] ]) + + return cls(B, field, **kwargs) + def inner_product(self, x, y): r""" Half of the trace inner product. @@ -1918,21 +2225,15 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra): paper:: sage: set_random_seed() - sage: n = ZZ.random_element(2,5) - sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular') - sage: B = M.transpose()*M - sage: J = BilinearFormEJA(n, B=B) + sage: J = BilinearFormEJA.random_instance() + sage: n = J.dimension() sage: x = J.random_element() sage: y = J.random_element() - sage: x.inner_product(y) == (x*y).trace()/2 + sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) True """ - xvec = x.to_vector() - xbar = xvec[1:] - yvec = y.to_vector() - ybar = yvec[1:] - return x[0]*y[0] + (self._B*xbar).inner_product(ybar) + return (self._B*x.to_vector()).inner_product(y.to_vector()) class JordanSpinEJA(BilinearFormEJA): @@ -1988,7 +2289,8 @@ class JordanSpinEJA(BilinearFormEJA): def __init__(self, n, field=AA, **kwargs): # This is a special case of the BilinearFormEJA with the identity # matrix as its bilinear form. - return super(JordanSpinEJA, self).__init__(n, field, **kwargs) + B = matrix.identity(field, n) + super(JordanSpinEJA, self).__init__(B, field, **kwargs) class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra): @@ -2022,8 +2324,198 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra): """ def __init__(self, field=AA, **kwargs): mult_table = [] - fdeja = super(TrivialEJA, self) + super(TrivialEJA, self).__init__(field, + mult_table, + check_axioms=False, + **kwargs) # The rank is zero using my definition, namely the dimension of the # largest subalgebra generated by any element. - fdeja.__init__(field, mult_table, **kwargs) self.rank.set_cache(0) + + @classmethod + def random_instance(cls, field=AA, **kwargs): + # We don't take a "size" argument so the superclass method is + # inappropriate for us. + return cls(field, **kwargs) + +class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): + 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. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: RealSymmetricEJA, + ....: DirectSumEJA) + + EXAMPLES:: + + sage: J1 = HadamardEJA(2) + sage: J2 = RealSymmetricEJA(3) + sage: J = DirectSumEJA(J1,J2) + sage: J.dimension() + 8 + sage: J.rank() + 5 + + """ + def __init__(self, J1, J2, field=AA, **kwargs): + 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(n) ] + for i in range(n) ] + for i in range(n1): + for j in range(n1): + p = (J1.monomial(i)*J1.monomial(j)).to_vector() + mult_table[i][j] = V(p.list() + [field.zero()]*n2) + + for i in range(n2): + for j in range(n2): + p = (J2.monomial(i)*J2.monomial(j)).to_vector() + mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list()) + + super(DirectSumEJA, self).__init__(field, + mult_table, + check_axioms=False, + **kwargs) + self.rank.set_cache(J1.rank() + J2.rank()) + + + def factors(self): + r""" + Return the pair of this algebra's factors. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: JordanSpinEJA, + ....: DirectSumEJA) + + EXAMPLES:: + + sage: J1 = HadamardEJA(2,QQ) + sage: J2 = JordanSpinEJA(3,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) + + """ + return self._factors + + def projections(self): + r""" + Return a pair of projections onto this algebra's factors. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + ....: ComplexHermitianEJA, + ....: DirectSumEJA) + + EXAMPLES:: + + 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: pi_left(J.one()).to_vector() + (1, 0) + sage: pi_right(J.one()).to_vector() + (1, 0, 0, 1) + + """ + (J1,J2) = self.factors() + n = J1.dimension() + pi_left = lambda x: J1.from_vector(x.to_vector()[:n]) + pi_right = lambda x: J2.from_vector(x.to_vector()[n:]) + return (pi_left, pi_right) + + def inclusions(self): + r""" + Return the pair of inclusion maps from our factors into us. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + ....: RealSymmetricEJA, + ....: DirectSumEJA) + + EXAMPLES:: + + sage: J1 = JordanSpinEJA(3) + sage: J2 = RealSymmetricEJA(2) + sage: J = DirectSumEJA(J1,J2) + sage: (iota_left, iota_right) = J.inclusions() + sage: iota_left(J1.zero()) == J.zero() + True + sage: iota_right(J2.zero()) == J.zero() + True + sage: J1.one().to_vector() + (1, 0, 0) + sage: iota_left(J1.one()).to_vector() + (1, 0, 0, 0, 0, 0) + sage: J2.one().to_vector() + (1, 0, 1) + sage: iota_right(J2.one()).to_vector() + (0, 0, 0, 1, 0, 1) + sage: J.one().to_vector() + (1, 0, 0, 1, 0, 1) + + """ + (J1,J2) = self.factors() + n = J1.dimension() + V_basis = self.vector_space().basis() + I1 = matrix.column(self.base_ring(), V_basis[:n]) + I2 = matrix.column(self.base_ring(), V_basis[n:]) + iota_left = lambda x: self.from_vector(I1*x.to_vector()) + iota_right = lambda x: self.from_vector(I2*+x.to_vector()) + return (iota_left, iota_right) + + 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. + + SETUP:: + + + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: QuaternionHermitianEJA, + ....: DirectSumEJA) + + EXAMPLE:: + + sage: J1 = HadamardEJA(3) + sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=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 + + """ + (pi_left, pi_right) = self.projections() + x1 = pi_left(x) + x2 = pi_right(x) + y1 = pi_left(y) + y2 = pi_right(y) + + return (x1.inner_product(y1) + x2.inner_product(y2))