X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=44d83147bbeceb4acae52d573749545dec3c4b99;hb=9708043704809263d2bd543de38c46b458b873cb;hp=3b5828fed55098dfd3f4c95bf2354ac4e38efc87;hpb=0fc6cf97abf8e091787ebae4e9cd60534ebdbc32;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 3b5828f..44d8314 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -24,15 +24,13 @@ 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 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): @@ -56,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 @@ -66,7 +64,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): mult_table, prefix='e', category=None, - natural_basis=None, + matrix_basis=None, check_field=True, check_axioms=True): """ @@ -117,7 +115,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): if not all( len(l) == n for l in mult_table ): raise ValueError("multiplication table is not square") - self._natural_basis = natural_basis + self._matrix_basis = matrix_basis if category is None: category = MagmaticAlgebras(field).FiniteDimensional() @@ -136,10 +134,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # long run to have the multiplication table be in terms of # algebra elements. We do this after calling the superclass # constructor so that from_vector() knows what to do. - self._multiplication_table = [ - list(map(lambda x: self.from_vector(x), ls)) - for ls in mult_table - ] + self._multiplication_table = [ [ self.vector_space().zero() + for i in range(n) ] + for j in range(n) ] + # take advantage of symmetry + for i in range(n): + for j in range(i+1): + elt = self.from_vector(mult_table[i][j]) + self._multiplication_table[i][j] = elt + self._multiplication_table[j][i] = elt if check_axioms: if not self._is_commutative(): @@ -151,7 +154,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): def _element_constructor_(self, elt): """ - Construct an element of this algebra from its natural + Construct an element of this algebra from its vector or matrix representation. This gets called only after the parent element _call_ method @@ -179,13 +182,13 @@ 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: Ensure that we can convert any element of the two non-matrix - simple algebras (whose natural representations are their usual - vector representations) back and forth faithfully:: + simple algebras (whose matrix representations are columns) + back and forth faithfully:: sage: set_random_seed() sage: J = HadamardEJA.random_instance() @@ -196,9 +199,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: x = J.random_element() sage: J(x.to_vector().column()) == x 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. @@ -210,20 +212,23 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # that the integer 3 belongs to the space of 2-by-2 matrices. raise ValueError(msg) - natural_basis = self.natural_basis() - basis_space = natural_basis[0].matrix_space() - if elt not in basis_space: + if elt not in self.matrix_space(): raise ValueError(msg) # Thanks for nothing! Matrix spaces aren't vector spaces in - # Sage, so we have to figure out its natural-basis coordinates + # Sage, so we have to figure out its matrix-basis coordinates # ourselves. We use the basis space's ring instead of the # element's ring because the basis space might be an algebraic # closure whereas the base ring of the 3-by-3 identity matrix # 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)) + V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols()) + W = V.span_of_basis( _mat2vec(s) for s in self.matrix_basis() ) + + 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): @@ -373,6 +378,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): """ @@ -384,7 +411,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: @@ -397,10 +426,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): @@ -464,18 +517,33 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return table(M, header_row=True, header_column=True, frame=True) - def natural_basis(self): + def matrix_basis(self): """ - Return a more-natural representation of this algebra's basis. + Return an (often more natural) representation of this algebras + basis as an ordered tuple of matrices. + + Every finite-dimensional Euclidean Jordan Algebra is a, up to + Jordan isomorphism, a direct sum of five simple + algebras---four of which comprise Hermitian matrices. And the + last type of algebra can of course be thought of as `n`-by-`1` + column matrices (ambiguusly called column vectors) to avoid + special cases. As a result, matrices (and column vectors) are + a natural representation format for Euclidean Jordan algebra + elements. - Every finite-dimensional Euclidean Jordan Algebra is a direct - sum of five simple algebras, four of which comprise Hermitian - matrices. This method returns the original "natural" basis - for our underlying vector space. (Typically, the natural basis - is used to construct the multiplication table in the first place.) + But, when we construct an algebra from a basis of matrices, + those matrix representations are lost in favor of coordinate + vectors *with respect to* that basis. We could eventually + convert back if we tried hard enough, but having the original + representations handy is valuable enough that we simply store + them and return them from this method. - Note that this will always return a matrix. The standard basis - in `R^n` will be returned as `n`-by-`1` column matrices. + Why implement this for non-matrix algebras? Avoiding special + cases for the :class:`BilinearFormEJA` pays with simplicity in + its own right. But mainly, we would like to be able to assume + that elements of a :class:`DirectSumEJA` can be displayed + nicely, without having to have special classes for direct sums + one of whose components was a matrix algebra. SETUP:: @@ -487,7 +555,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = RealSymmetricEJA(2) sage: J.basis() Finite family {0: e0, 1: e1, 2: e2} - sage: J.natural_basis() + sage: J.matrix_basis() ( [1 0] [ 0 0.7071067811865475?] [0 0] [0 0], [0.7071067811865475? 0], [0 1] @@ -498,50 +566,38 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = JordanSpinEJA(2) sage: J.basis() Finite family {0: e0, 1: e1} - sage: J.natural_basis() + sage: J.matrix_basis() ( [1] [0] [0], [1] ) - """ - if self._natural_basis is None: - M = self.natural_basis_space() + if self._matrix_basis is None: + M = self.matrix_space() return tuple( M(b.to_vector()) for b in self.basis() ) else: - return self._natural_basis + return self._matrix_basis - def natural_basis_space(self): + def matrix_space(self): """ - Return the matrix space in which this algebra's natural basis - elements live. + Return the matrix space in which this algebra's elements live, if + we think of them as matrices (including column vectors of the + appropriate size). 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. + (where `n` is zero), to ensure that two elements of the matrix + space (empty matrices) can be multiplied. + + Matrix algebras override this with something more useful. """ if self.is_trivial(): return MatrixSpace(self.base_ring(), 0) - elif self._natural_basis is None or len(self._natural_basis) == 0: + elif self._matrix_basis is None or len(self._matrix_basis) == 0: return MatrixSpace(self.base_ring(), self.dimension(), 1) else: - 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() + return self._matrix_basis[0].matrix_space() @cached_method @@ -668,22 +724,22 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Vector space of degree 6 and dimension 2... sage: J1 Euclidean Jordan algebra of dimension 3... - sage: J0.one().natural_representation() + sage: J0.one().to_matrix() [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() + sage: J.from_vector(J5.basis()[0]).to_matrix() [ 0 0 1/2*sqrt(2)] [ 0 0 0] [1/2*sqrt(2) 0 0] - sage: J.from_vector(J5.basis()[1]).natural_representation() + sage: J.from_vector(J5.basis()[1]).to_matrix() [ 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() + sage: J1.one().to_matrix() [1 0 0] [0 1 0] [0 0 0] @@ -733,6 +789,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 @@ -841,8 +899,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() @@ -937,38 +994,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Ensure that computing the rank actually works, since the ranks of all simple algebras are known and will be cached by default:: - sage: J = HadamardEJA(4) - sage: J.rank.clear_cache() - sage: J.rank() - 4 - - :: - - sage: J = JordanSpinEJA(4) - sage: J.rank.clear_cache() - sage: J.rank() - 2 - - :: - - sage: J = RealSymmetricEJA(3) - sage: J.rank.clear_cache() - sage: J.rank() - 3 - - :: - - sage: J = ComplexHermitianEJA(2) - sage: J.rank.clear_cache() - sage: J.rank() - 2 - - :: + sage: set_random_seed() # long time + sage: J = random_eja() # long time + sage: caches = J.rank() # long time + sage: J.rank.clear_cache() # long time + sage: J.rank() == cached # long time + True - sage: J = QuaternionHermitianEJA(2) - sage: J.rank.clear_cache() - sage: J.rank() - 2 """ return len(self._charpoly_coefficients()) @@ -993,6 +1025,75 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Element = FiniteDimensionalEuclideanJordanAlgebraElement +class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra): + def __init__(self, + field, + basis, + jordan_product, + inner_product, + orthonormalize=True, + prefix='e', + category=None, + check_field=True, + check_axioms=True): + + n = len(basis) + vector_basis = basis + + from sage.matrix.matrix import is_Matrix + basis_is_matrices = False + + degree = 0 + if n > 0: + if is_Matrix(basis[0]): + basis_is_matrices = True + vector_basis = tuple( map(_mat2vec,basis) ) + degree = basis[0].nrows()**2 + else: + degree = basis[0].degree() + + V = VectorSpace(field, degree) + + self._deorthonormalization_matrix = matrix.identity(field,n) + if orthonormalize: + A = matrix(field, vector_basis) + # uh oh, this is only the "usual" inner product + Q,R = A.gram_schmidt(orthonormal=True) + self._deorthonormalization_matrix = R.inverse().transpose() + vector_basis = Q.rows() + W = V.span_of_basis( vector_basis ) + if basis_is_matrices: + from mjo.eja.eja_utils import _vec2mat + basis = tuple( map(_vec2mat,vector_basis) ) + + mult_table = [ [0 for i in range(n)] for j in range(n) ] + ip_table = [ [0 for i in range(n)] for j in range(n) ] + + for i in range(n): + for j in range(i+1): + # do another mat2vec because the multiplication + # table is in terms of vectors + elt = _mat2vec(jordan_product(basis[i],basis[j])) + elt = W.coordinate_vector(elt) + mult_table[i][j] = elt + mult_table[j][i] = elt + ip = inner_product(basis[i],basis[j]) + ip_table[i][j] = ip + ip_table[j][i] = ip + + self._inner_product_matrix = matrix(field,ip_table) + + if basis_is_matrices: + for m in basis: + m.set_immutable() + + super().__init__(field, + mult_table, + prefix, + category, + basis, + check_field, + check_axioms) class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): r""" @@ -1038,7 +1139,7 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr return superclass._charpoly_coefficients() mult_table = tuple( - map(lambda x: x.to_vector(), ls) + tuple(map(lambda x: x.to_vector(), ls)) for ls in self._multiplication_table ) @@ -1068,26 +1169,25 @@ class ConcreteEuclideanJordanAlgebra: TESTS: - Our natural basis is normalized with respect to the natural inner - product unless we specify otherwise:: + Our basis is normalized with respect to the algebra's inner + product, unless we specify otherwise:: sage: set_random_seed() sage: J = ConcreteEuclideanJordanAlgebra.random_instance() sage: all( b.norm() == 1 for b in J.gens() ) True - Since our natural basis is normalized with respect to the natural - inner product, and since we know that this algebra is an EJA, any + Since our basis is orthonormal with respect to the algebra's inner + product, and since we know that this algebra is an EJA, any left-multiplication operator's matrix will be symmetric because - natural->EJA basis representation is an isometry and within the EJA - the operator is self-adjoint by the Jordan axiom:: + natural->EJA basis representation is an isometry and within the + EJA the operator is self-adjoint by the Jordan axiom:: 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 - """ @staticmethod @@ -1134,6 +1234,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 @@ -1145,14 +1249,41 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt()) basis = tuple( s.change_ring(field) for s in basis ) self._basis_normalizers = tuple( - ~(self.natural_inner_product(s,s).sqrt()) for s in basis ) + ~(self.matrix_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.matrix_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, - natural_basis=basis, + mult_table, + matrix_basis=basis, **kwargs) if algebra_dim == 0: @@ -1175,7 +1306,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): # entries in a nice field already. Just compute the thing. return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients() - basis = ( (b/n) for (b,n) in zip(self.natural_basis(), + basis = ( (b/n) for (b,n) in zip(self.matrix_basis(), self._basis_normalizers) ) # Do this over the rationals and convert back at the end. @@ -1210,39 +1341,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): """ @@ -1267,9 +1365,8 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): """ raise NotImplementedError - @classmethod - def natural_inner_product(cls,X,Y): + def matrix_inner_product(cls,X,Y): Xu = cls.real_unembed(X) Yu = cls.real_unembed(Y) tr = (Xu*Yu).trace() @@ -1348,9 +1445,9 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, sage: set_random_seed() sage: J = RealSymmetricEJA.random_instance() sage: x,y = J.random_elements(2) - sage: actual = (x*y).natural_representation() - sage: X = x.natural_representation() - sage: Y = y.natural_representation() + sage: actual = (x*y).to_matrix() + sage: X = x.to_matrix() + sage: Y = y.to_matrix() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True @@ -1549,9 +1646,9 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @classmethod - def natural_inner_product(cls,X,Y): + def matrix_inner_product(cls,X,Y): """ - Compute a natural inner product in this algebra directly from + Compute a matrix inner product in this algebra directly from its real embedding. SETUP:: @@ -1566,17 +1663,17 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() sage: x,y = J.random_elements(2) - sage: Xe = x.natural_representation() - sage: Ye = y.natural_representation() + sage: Xe = x.to_matrix() + sage: Ye = y.to_matrix() sage: X = ComplexHermitianEJA.real_unembed(Xe) sage: Y = ComplexHermitianEJA.real_unembed(Ye) sage: expected = (X*Y).trace().real() - sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye) + sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye) sage: actual == expected True """ - return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2 + return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/2 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, @@ -1617,9 +1714,9 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() sage: x,y = J.random_elements(2) - sage: actual = (x*y).natural_representation() - sage: X = x.natural_representation() - sage: Y = y.natural_representation() + sage: actual = (x*y).to_matrix() + sage: X = x.to_matrix() + sage: Y = y.to_matrix() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True @@ -1844,9 +1941,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @classmethod - def natural_inner_product(cls,X,Y): + def matrix_inner_product(cls,X,Y): """ - Compute a natural inner product in this algebra directly from + Compute a matrix inner product in this algebra directly from its real embedding. SETUP:: @@ -1861,17 +1958,17 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): sage: set_random_seed() sage: J = QuaternionHermitianEJA.random_instance() sage: x,y = J.random_elements(2) - sage: Xe = x.natural_representation() - sage: Ye = y.natural_representation() + sage: Xe = x.to_matrix() + sage: Ye = y.to_matrix() sage: X = QuaternionHermitianEJA.real_unembed(Xe) sage: Y = QuaternionHermitianEJA.real_unembed(Ye) sage: expected = (X*Y).trace().coefficient_tuple()[0] - sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye) + sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye) sage: actual == expected True """ - return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4 + return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/4 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, @@ -1912,9 +2009,9 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, sage: set_random_seed() sage: J = QuaternionHermitianEJA.random_instance() sage: x,y = J.random_elements(2) - sage: actual = (x*y).natural_representation() - sage: X = x.natural_representation() - sage: Y = y.natural_representation() + sage: actual = (x*y).to_matrix() + sage: X = x.to_matrix() + sage: Y = y.to_matrix() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True @@ -2058,6 +2155,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, @@ -2085,31 +2192,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""" @@ -2191,7 +2273,6 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra, True """ def __init__(self, B, field=AA, **kwargs): - self._B = B n = B.nrows() if not B.is_positive_definite(): @@ -2212,13 +2293,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: @@ -2257,35 +2349,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): """ @@ -2331,9 +2394,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 """ @@ -2393,6 +2456,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, @@ -2418,7 +2482,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): SETUP:: - sage: from mjo.eja.eja_algebra import (HadamardEJA, + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, ....: RealSymmetricEJA, ....: DirectSumEJA) @@ -2432,8 +2497,25 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: J.rank() 5 + TESTS: + + The external direct sum construction is only valid when the two factors + have the same base ring; an error is raised otherwise:: + + sage: set_random_seed() + sage: J1 = random_eja(AA) + sage: J2 = random_eja(QQ) + sage: J = DirectSumEJA(J1,J2) + Traceback (most recent call last): + ... + ValueError: algebras must share the same base field + """ - def __init__(self, J1, J2, field=AA, **kwargs): + def __init__(self, J1, J2, **kwargs): + if J1.base_ring() != J2.base_ring(): + raise ValueError("algebras must share the same base field") + field = J1.base_ring() + self._factors = (J1, J2) n1 = J1.dimension() n2 = J2.dimension() @@ -2505,9 +2587,16 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): """ (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:]) + m = J1.dimension() + n = J2.dimension() + V_basis = self.vector_space().basis() + # Need to specify the dimensions explicitly so that we don't + # wind up with a zero-by-zero matrix when we want e.g. a + # zero-by-two matrix (important for composing things). + P1 = matrix(self.base_ring(), m, m+n, V_basis[:m]) + P2 = matrix(self.base_ring(), n, m+n, V_basis[m:]) + pi_left = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J1,P1) + pi_right = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J2,P2) return (pi_left, pi_right) def inclusions(self): @@ -2516,7 +2605,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + sage: from mjo.eja.eja_algebra import (random_eja, + ....: JordanSpinEJA, ....: RealSymmetricEJA, ....: DirectSumEJA) @@ -2541,14 +2631,39 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: J.one().to_vector() (1, 0, 0, 1, 0, 1) + TESTS: + + Composing a projection with the corresponding inclusion should + produce the identity map, and mismatching them should produce + the zero map:: + + sage: set_random_seed() + sage: J1 = random_eja() + sage: J2 = random_eja() + sage: J = DirectSumEJA(J1,J2) + sage: (iota_left, iota_right) = J.inclusions() + sage: (pi_left, pi_right) = J.projections() + sage: pi_left*iota_left == J1.one().operator() + True + sage: pi_right*iota_right == J2.one().operator() + True + sage: (pi_left*iota_right).is_zero() + True + sage: (pi_right*iota_left).is_zero() + True + """ (J1,J2) = self.factors() - n = J1.dimension() + m = J1.dimension() + n = J2.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()) + # Need to specify the dimensions explicitly so that we don't + # wind up with a zero-by-zero matrix when we want e.g. a + # two-by-zero matrix (important for composing things). + I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m]) + I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:]) + iota_left = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,self,I1) + iota_right = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,self,I2) return (iota_left, iota_right) def inner_product(self, x, y): @@ -2567,7 +2682,7 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): EXAMPLE:: - sage: J1 = HadamardEJA(3) + sage: J1 = HadamardEJA(3,QQ) sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False) sage: J = DirectSumEJA(J1,J2) sage: x1 = J1.one()