X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=2ea1a9267623d95d02231654a48a470efc186a40;hb=4c6ab92d378613a9a053b62382b7e0cda7c450ab;hp=6252cc29a729472ef7f182cd3393ab3db98b6b58;hpb=39d8d3190b721ea21e0e86618d774437bc1eeb35;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 6252cc2..2ea1a92 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -24,15 +24,12 @@ from sage.combinat.free_module import CombinatorialFreeModule from sage.matrix.constructor import matrix from sage.matrix.matrix_space import MatrixSpace from sage.misc.cachefunc import cached_method -from sage.misc.lazy_import import lazy_import from sage.misc.table import table from sage.modules.free_module import FreeModule, VectorSpace from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF, PolynomialRing, QuadraticField) from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement -lazy_import('mjo.eja.eja_subalgebra', - 'FiniteDimensionalEuclideanJordanSubalgebra') from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator from mjo.eja.eja_utils import _mat2vec @@ -57,7 +54,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J(1) Traceback (most recent call last): ... - ValueError: not a naturally-represented algebra element + ValueError: not an element of this algebra """ return None @@ -180,7 +177,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J(A) Traceback (most recent call last): ... - ArithmeticError: vector is not in free module + ValueError: not an element of this algebra TESTS: @@ -199,7 +196,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): True """ - msg = "not a naturally-represented algebra element" + msg = "not an element of this algebra" if elt == 0: # The superclass implementation of random_element() # needs to be able to coerce "0" into the algebra. @@ -224,7 +221,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # could be QQ instead of QQbar. V = VectorSpace(basis_space.base_ring(), elt.nrows()*elt.ncols()) W = V.span_of_basis( _mat2vec(s) for s in natural_basis ) - coords = W.coordinate_vector(_mat2vec(elt)) + + try: + coords = W.coordinate_vector(_mat2vec(elt)) + except ArithmeticError: # vector is not in free module + raise ValueError(msg) + return self.from_vector(coords) def _repr_(self): @@ -374,6 +376,28 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return (t**r + sum( a[k]*(t**k) for k in range(r) )) + def coordinate_polynomial_ring(self): + r""" + The multivariate polynomial ring in which this algebra's + :meth:`characteristic_polynomial_of` lives. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: RealSymmetricEJA) + + EXAMPLES:: + + sage: J = HadamardEJA(2) + sage: J.coordinate_polynomial_ring() + Multivariate Polynomial Ring in X1, X2... + sage: J = RealSymmetricEJA(3,QQ) + sage: J.coordinate_polynomial_ring() + Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6... + + """ + var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) ) + return PolynomialRing(self.base_ring(), var_names) def inner_product(self, x, y): """ @@ -385,7 +409,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): SETUP:: - sage: from mjo.eja.eja_algebra import random_eja + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, + ....: BilinearFormEJA) EXAMPLES: @@ -398,10 +424,34 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: (x*y).inner_product(z) == y.inner_product(x*z) True + 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: actual = x.inner_product(y) + sage: expected = x.to_vector().inner_product(y.to_vector()) + sage: actual == expected + True + + Ensure that this is one-half of the trace inner-product in a + BilinearFormEJA that isn't just the reals (when ``n`` isn't + one). This is in Faraut and Koranyi, and also my "On the + symmetry..." paper:: + + sage: set_random_seed() + sage: J = BilinearFormEJA.random_instance() + sage: n = J.dimension() + sage: x = J.random_element() + sage: y = J.random_element() + sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) + True """ - X = x.natural_representation() - Y = y.natural_representation() - return self.natural_inner_product(X,Y) + B = self._inner_product_matrix + return (B*x.to_vector()).inner_product(y.to_vector()) def is_trivial(self): @@ -531,20 +581,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return self._natural_basis[0].matrix_space() - @staticmethod - def natural_inner_product(X,Y): - """ - Compute the inner product of two naturally-represented elements. - - For example in the real symmetric matrix EJA, this will compute - the trace inner-product of two n-by-n symmetric matrices. The - default should work for the real cartesian product EJA, the - Jordan spin EJA, and the real symmetric matrices. The others - will have to be overridden. - """ - return (X.conjugate_transpose()*Y).trace() - - @cached_method def one(self): """ @@ -734,6 +770,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): if not c.is_idempotent(): raise ValueError("element is not idempotent: %s" % c) + from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra + # Default these to what they should be if they turn out to be # trivial, because eigenspaces_left() won't return eigenvalues # corresponding to trivial spaces (e.g. it returns only the @@ -842,8 +880,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): of" function. """ n = self.dimension() - var_names = [ "X" + str(z) for z in range(1,n+1) ] - R = PolynomialRing(self.base_ring(), var_names) + R = self.coordinate_polynomial_ring() vars = R.gens() F = R.fraction_field() @@ -1086,7 +1123,7 @@ class ConcreteEuclideanJordanAlgebra: sage: set_random_seed() sage: J = ConcreteEuclideanJordanAlgebra.random_instance() sage: x = J.random_element() - sage: x.operator().matrix().is_symmetric() + sage: x.operator().is_self_adjoint() True """ @@ -1135,6 +1172,10 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): basis = tuple(basis) algebra_dim = len(basis) + degree = 0 # size of the matrices + if algebra_dim > 0: + degree = basis[0].nrows() + if algebra_dim > 1 and normalize_basis: # We'll need sqrt(2) to normalize the basis, and this # winds up in the multiplication table, so the whole @@ -1149,10 +1190,37 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): ~(self.natural_inner_product(s,s).sqrt()) for s in basis ) basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers)) - Qs = self.multiplication_table_from_matrix_basis(basis) + # Now compute the multiplication and inner product tables. + # We have to do this *after* normalizing the basis, because + # scaling affects the answers. + V = VectorSpace(field, degree**2) + W = V.span_of_basis( _mat2vec(s) for s in basis ) + mult_table = [[W.zero() for j in range(algebra_dim)] + for i in range(algebra_dim)] + ip_table = [[W.zero() for j in range(algebra_dim)] + for i in range(algebra_dim)] + for i in range(algebra_dim): + for j in range(algebra_dim): + mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2 + mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry)) + + try: + # HACK: ignore the error here if we don't need the + # inner product (as is the case when we construct + # a dummy QQ-algebra for fast charpoly coefficients. + ip_table[i][j] = self.natural_inner_product(basis[i], + basis[j]) + except: + pass + + try: + # HACK PART DEUX + self._inner_product_matrix = matrix(field,ip_table) + except: + pass super(MatrixEuclideanJordanAlgebra, self).__init__(field, - Qs, + mult_table, natural_basis=basis, **kwargs) @@ -1211,39 +1279,6 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): return tuple( a_i.subs(subs_dict) for a_i in a ) - @staticmethod - def multiplication_table_from_matrix_basis(basis): - """ - At least three of the five simple Euclidean Jordan algebras have the - symmetric multiplication (A,B) |-> (AB + BA)/2, where the - multiplication on the right is matrix multiplication. Given a basis - for the underlying matrix space, this function returns a - multiplication table (obtained by looping through the basis - elements) for an algebra of those matrices. - """ - # In S^2, for example, we nominally have four coordinates even - # though the space is of dimension three only. The vector space V - # 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() - - V = VectorSpace(field, dimension**2) - W = V.span_of_basis( _mat2vec(s) for s in basis ) - n = len(basis) - mult_table = [[W.zero() for j in range(n)] for i in range(n)] - for i in range(n): - for j in range(n): - mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2 - mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry)) - - return mult_table - - @staticmethod def real_embed(M): """ @@ -1268,7 +1303,6 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): """ raise NotImplementedError - @classmethod def natural_inner_product(cls,X,Y): Xu = cls.real_unembed(X) @@ -2059,6 +2093,16 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra, mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ] for i in range(n) ] + # Inner products are real numbers and not algebra + # elements, so once we turn the algebra element + # into a vector in inner_product(), we never go + # back. As a result -- contrary to what we do with + # self._multiplication_table -- we store the inner + # product table as a plain old matrix and not as + # an algebra operator. + ip_table = matrix.identity(field,n) + self._inner_product_matrix = ip_table + super(HadamardEJA, self).__init__(field, mult_table, check_axioms=False, @@ -2086,31 +2130,6 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra, return cls(n, field, **kwargs) - 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, ConcreteEuclideanJordanAlgebra): r""" @@ -2192,7 +2211,6 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra, True """ def __init__(self, B, field=AA, **kwargs): - self._B = B n = B.nrows() if not B.is_positive_definite(): @@ -2213,13 +2231,24 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra, z = V([z0] + zbar.list()) mult_table[i][j] = z - # 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). + # Inner products are real numbers and not algebra + # elements, so once we turn the algebra element + # into a vector in inner_product(), we never go + # back. As a result -- contrary to what we do with + # self._multiplication_table -- we store the inner + # product table as a plain old matrix and not as + # an algebra operator. + ip_table = B + self._inner_product_matrix = ip_table + super(BilinearFormEJA, self).__init__(field, mult_table, check_axioms=False, **kwargs) + + # 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). self.rank.set_cache(min(n,2)) if n == 0: @@ -2258,35 +2287,6 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra, return cls(B, field, **kwargs) - def inner_product(self, x, y): - r""" - Half of the trace inner product. - - This is defined so that the special case of the Jordan spin - algebra gets the usual inner product. - - SETUP:: - - sage: from mjo.eja.eja_algebra import BilinearFormEJA - - TESTS: - - Ensure that this is one-half of the trace inner-product when - the algebra isn't just the reals (when ``n`` isn't one). This - is in Faraut and Koranyi, and also my "On the symmetry..." - paper:: - - sage: set_random_seed() - sage: J = BilinearFormEJA.random_instance() - sage: n = J.dimension() - sage: x = J.random_element() - sage: y = J.random_element() - sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) - True - - """ - return (self._B*x.to_vector()).inner_product(y.to_vector()) - class JordanSpinEJA(BilinearFormEJA): """ @@ -2332,9 +2332,9 @@ class JordanSpinEJA(BilinearFormEJA): sage: set_random_seed() sage: J = JordanSpinEJA.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) + sage: actual = x.inner_product(y) + sage: expected = x.to_vector().inner_product(y.to_vector()) + sage: actual == expected True """ @@ -2394,6 +2394,7 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, """ def __init__(self, field=AA, **kwargs): mult_table = [] + self._inner_product_matrix = matrix(field,0) super(TrivialEJA, self).__init__(field, mult_table, check_axioms=False,