X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feuclidean_jordan_algebra.py;h=b25839815289ba6f60630a14b790dbfe29f8969b;hb=5978697af0b924d2d2b6174975bf06f1830746f3;hp=c0b7787a535b7754c446fcb0046c300b21695579;hpb=6ca4fb6c40e9e8f65eec8cafd0044992fae79f02;p=sage.d.git diff --git a/mjo/eja/euclidean_jordan_algebra.py b/mjo/eja/euclidean_jordan_algebra.py index c0b7787..b258398 100644 --- a/mjo/eja/euclidean_jordan_algebra.py +++ b/mjo/eja/euclidean_jordan_algebra.py @@ -20,7 +20,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): names='e', assume_associative=False, category=None, - rank=None): + rank=None, + natural_basis=None, + inner_product=None): n = len(mult_table) mult_table = [b.base_extend(field) for b in mult_table] for b in mult_table: @@ -43,7 +45,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): assume_associative=assume_associative, names=names, category=cat, - rank=rank) + rank=rank, + natural_basis=natural_basis, + inner_product=inner_product) def __init__(self, field, @@ -51,7 +55,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): names='e', assume_associative=False, category=None, - rank=None): + rank=None, + natural_basis=None, + inner_product=None): """ EXAMPLES: @@ -66,6 +72,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): """ self._rank = rank + self._natural_basis = natural_basis + self._inner_product = inner_product fda = super(FiniteDimensionalEuclideanJordanAlgebra, self) fda.__init__(field, mult_table, @@ -80,6 +88,63 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): fmt = "Euclidean Jordan algebra of degree {} over {}" return fmt.format(self.degree(), self.base_ring()) + + def inner_product(self, x, y): + """ + The inner product associated with this Euclidean Jordan algebra. + + Will default to the trace inner product if nothing else. + """ + if (not x in self) or (not y in self): + raise TypeError("arguments must live in this algebra") + if self._inner_product is None: + return x.trace_inner_product(y) + else: + return self._inner_product(x,y) + + + def natural_basis(self): + """ + Return a more-natural representation of this algebra's basis. + + 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.) + + Note that this will always return a matrix. The standard basis + in `R^n` will be returned as `n`-by-`1` column matrices. + + EXAMPLES:: + + sage: J = RealSymmetricSimpleEJA(2) + sage: J.basis() + Family (e0, e1, e2) + sage: J.natural_basis() + ( + [1 0] [0 1] [0 0] + [0 0], [1 0], [0 1] + ) + + :: + + sage: J = JordanSpinSimpleEJA(2) + sage: J.basis() + Family (e0, e1) + sage: J.natural_basis() + ( + [1] [0] + [0], [1] + ) + + """ + if self._natural_basis is None: + return tuple( b.vector().column() for b in self.basis() ) + else: + return self._natural_basis + + def rank(self): """ Return the rank of this EJA. @@ -95,6 +160,51 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): An element of a Euclidean Jordan algebra. """ + def __init__(self, A, elt=None): + """ + EXAMPLES: + + The identity in `S^n` is converted to the identity in the EJA:: + + sage: J = RealSymmetricSimpleEJA(3) + sage: I = identity_matrix(QQ,3) + sage: J(I) == J.one() + True + + This skew-symmetric matrix can't be represented in the EJA:: + + sage: J = RealSymmetricSimpleEJA(3) + sage: A = matrix(QQ,3, lambda i,j: i-j) + sage: J(A) + Traceback (most recent call last): + ... + ArithmeticError: vector is not in free module + + """ + # Goal: if we're given a matrix, and if it lives in our + # parent algebra's "natural ambient space," convert it + # into an algebra element. + # + # The catch is, we make a recursive call after converting + # the given matrix into a vector that lives in the algebra. + # This we need to try the parent class initializer first, + # to avoid recursing forever if we're given something that + # already fits into the algebra, but also happens to live + # in the parent's "natural ambient space" (this happens with + # vectors in R^n). + try: + FiniteDimensionalAlgebraElement.__init__(self, A, elt) + except ValueError: + natural_basis = A.natural_basis() + if elt in natural_basis[0].matrix_space(): + # Thanks for nothing! Matrix spaces aren't vector + # spaces in Sage, so we have to figure out its + # natural-basis coordinates ourselves. + V = VectorSpace(elt.base_ring(), elt.nrows()**2) + W = V.span( _mat2vec(s) for s in natural_basis ) + coords = W.coordinates(_mat2vec(elt)) + FiniteDimensionalAlgebraElement.__init__(self, A, coords) + def __pow__(self, n): """ Return ``self`` raised to the power ``n``. @@ -161,6 +271,62 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): raise NotImplementedError('irregular element') + def inner_product(self, other): + """ + Return the parent algebra's inner product of myself and ``other``. + + EXAMPLES: + + The inner product in the Jordan spin algebra is the usual + inner product on `R^n` (this example only works because the + basis for the Jordan algebra is the standard basis in `R^n`):: + + sage: J = JordanSpinSimpleEJA(3) + sage: x = vector(QQ,[1,2,3]) + sage: y = vector(QQ,[4,5,6]) + sage: x.inner_product(y) + 32 + sage: J(x).inner_product(J(y)) + 32 + + The inner product on `S^n` is ` = trace(X*Y)`, where + multiplication is the usual matrix multiplication in `S^n`, + so the inner product of the identity matrix with itself + should be the `n`:: + + sage: J = RealSymmetricSimpleEJA(3) + sage: J.one().inner_product(J.one()) + 3 + + Likewise, the inner product on `C^n` is ` = + Re(trace(X*Y))`, where we must necessarily take the real + part because the product of Hermitian matrices may not be + Hermitian:: + + sage: J = ComplexHermitianSimpleEJA(3) + sage: J.one().inner_product(J.one()) + 3 + + TESTS: + + Ensure that we can always compute an inner product, and that + it gives us back a real number:: + + sage: set_random_seed() + sage: J = random_eja() + sage: x = J.random_element() + sage: y = J.random_element() + sage: x.inner_product(y) in RR + True + + """ + P = self.parent() + if not other in P: + raise TypeError("'other' must live in the same algebra") + + return P.inner_product(self, other) + + def operator_commutes_with(self, other): """ Return whether or not this element operator-commutes @@ -191,7 +357,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): """ if not other in self.parent(): - raise ArgumentError("'other' must live in the same algebra") + raise TypeError("'other' must live in the same algebra") A = self.operator_matrix() B = other.operator_matrix() @@ -281,7 +447,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): # TODO: we can do better once the call to is_invertible() # doesn't crash on irregular elements. #if not self.is_invertible(): - # raise ArgumentError('element is not invertible') + # raise ValueError('element is not invertible') # We do this a little different than the usual recursive # call to a finite-dimensional algebra element, because we @@ -414,7 +580,93 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): return self.span_of_powers().dimension() - def matrix(self): + def minimal_polynomial(self): + """ + EXAMPLES:: + + sage: set_random_seed() + sage: x = random_eja().random_element() + sage: x.degree() == x.minimal_polynomial().degree() + True + + :: + + sage: set_random_seed() + sage: x = random_eja().random_element() + sage: x.degree() == x.minimal_polynomial().degree() + True + + The minimal polynomial and the characteristic polynomial coincide + and are known (see Alizadeh, Example 11.11) for all elements of + the spin factor algebra that aren't scalar multiples of the + identity:: + + sage: set_random_seed() + sage: n = ZZ.random_element(2,10) + sage: J = JordanSpinSimpleEJA(n) + sage: y = J.random_element() + sage: while y == y.coefficient(0)*J.one(): + ....: y = J.random_element() + sage: y0 = y.vector()[0] + sage: y_bar = y.vector()[1:] + sage: actual = y.minimal_polynomial() + sage: x = SR.symbol('x', domain='real') + sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2) + sage: bool(actual == expected) + True + + """ + # The element we're going to call "minimal_polynomial()" on. + # Either myself, interpreted as an element of a finite- + # dimensional algebra, or an element of an associative + # subalgebra. + elt = None + + if self.parent().is_associative(): + elt = FiniteDimensionalAlgebraElement(self.parent(), self) + else: + V = self.span_of_powers() + assoc_subalg = self.subalgebra_generated_by() + # Mis-design warning: the basis used for span_of_powers() + # and subalgebra_generated_by() must be the same, and in + # the same order! + elt = assoc_subalg(V.coordinates(self.vector())) + + # Recursive call, but should work since elt lives in an + # associative algebra. + return elt.minimal_polynomial() + + + def natural_representation(self): + """ + Return a more-natural representation of this element. + + 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" representation of this element as a Hermitian + matrix, if it has one. If not, you get the usual representation. + + EXAMPLES:: + + sage: J = ComplexHermitianSimpleEJA(3) + sage: J.one() + e0 + e5 + e8 + sage: J.one().natural_representation() + [1 0 0 0 0 0] + [0 1 0 0 0 0] + [0 0 1 0 0 0] + [0 0 0 1 0 0] + [0 0 0 0 1 0] + [0 0 0 0 0 1] + + """ + B = self.parent().natural_basis() + W = B[0].matrix_space() + return W.linear_combination(zip(self.vector(), B)) + + + def operator_matrix(self): """ Return the matrix that represents left- (or right-) multiplication by this element in the parent algebra. @@ -480,71 +732,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self) return fda_elt.matrix().transpose() - # - # The plan is to eventually phase out "matrix()", which sounds - # too much like "matrix_representation()", in favor of the more- - # accurate "operator_matrix()". But we need to override matrix() - # to keep parent class methods happy in the meantime. - # - operator_matrix = matrix - - - def minimal_polynomial(self): - """ - EXAMPLES:: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: x.degree() == x.minimal_polynomial().degree() - True - - :: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: x.degree() == x.minimal_polynomial().degree() - True - - The minimal polynomial and the characteristic polynomial coincide - and are known (see Alizadeh, Example 11.11) for all elements of - the spin factor algebra that aren't scalar multiples of the - identity:: - - sage: set_random_seed() - sage: n = ZZ.random_element(2,10) - sage: J = JordanSpinSimpleEJA(n) - sage: y = J.random_element() - sage: while y == y.coefficient(0)*J.one(): - ....: y = J.random_element() - sage: y0 = y.vector()[0] - sage: y_bar = y.vector()[1:] - sage: actual = y.minimal_polynomial() - sage: x = SR.symbol('x', domain='real') - sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2) - sage: bool(actual == expected) - True - - """ - # The element we're going to call "minimal_polynomial()" on. - # Either myself, interpreted as an element of a finite- - # dimensional algebra, or an element of an associative - # subalgebra. - elt = None - - if self.parent().is_associative(): - elt = FiniteDimensionalAlgebraElement(self.parent(), self) - else: - V = self.span_of_powers() - assoc_subalg = self.subalgebra_generated_by() - # Mis-design warning: the basis used for span_of_powers() - # and subalgebra_generated_by() must be the same, and in - # the same order! - elt = assoc_subalg(V.coordinates(self.vector())) - - # Recursive call, but should work since elt lives in an - # associative algebra. - return elt.minimal_polynomial() - def quadratic_representation(self, other=None): """ @@ -616,7 +803,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): if other is None: other=self elif not other in self.parent(): - raise ArgumentError("'other' must live in the same algebra") + raise TypeError("'other' must live in the same algebra") L = self.operator_matrix() M = other.operator_matrix() @@ -779,7 +966,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): Return the trace inner product of myself and ``other``. """ if not other in self.parent(): - raise ArgumentError("'other' must live in the same algebra") + raise TypeError("'other' must live in the same algebra") return (self*other).trace() @@ -817,7 +1004,10 @@ def eja_rn(dimension, field=QQ): Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i)) for i in xrange(dimension) ] - return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension) + return FiniteDimensionalEuclideanJordanAlgebra(field, + Qs, + rank=dimension, + inner_product=_usual_ip) @@ -872,7 +1062,7 @@ def _real_symmetric_basis(n, field=QQ): # Beware, orthogonal but not normalized! Sij = Eij + Eij.transpose() S.append(Sij) - return S + return tuple(S) def _complex_hermitian_basis(n, field=QQ): @@ -909,9 +1099,15 @@ def _complex_hermitian_basis(n, field=QQ): S.append(Sij_real) Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose()) S.append(Sij_imag) - return S + return tuple(S) +def _mat2vec(m): + return vector(m.base_ring(), m.list()) + +def _vec2mat(v): + return matrix(v.base_ring(), sqrt(v.degree()), v.list()) + def _multiplication_table_from_matrix_basis(basis): """ At least three of the five simple Euclidean Jordan algebras have the @@ -919,7 +1115,10 @@ def _multiplication_table_from_matrix_basis(basis): 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. + elements) for an algebra of those matrices. A reordered copy + of the basis is also returned to work around the fact that + the ``span()`` in this function will change the order of the basis + from what we think it is, to... something else. """ # In S^2, for example, we nominally have four coordinates even # though the space is of dimension three only. The vector space V @@ -929,19 +1128,13 @@ def _multiplication_table_from_matrix_basis(basis): field = basis[0].base_ring() dimension = basis[0].nrows() - def mat2vec(m): - return vector(field, m.list()) - - def vec2mat(v): - return matrix(field, dimension, v.list()) - V = VectorSpace(field, dimension**2) - W = V.span( mat2vec(s) for s in basis ) + W = V.span( _mat2vec(s) for s in basis ) # Taking the span above reorders our basis (thanks, jerk!) so we # need to put our "matrix basis" in the same order as the # (reordered) vector basis. - S = [ vec2mat(b) for b in W.basis() ] + S = tuple( _vec2mat(b) for b in W.basis() ) Qs = [] for s in S: @@ -954,12 +1147,12 @@ def _multiplication_table_from_matrix_basis(basis): # why we're computing rows here and not columns. Q_rows = [] for t in S: - this_row = mat2vec((s*t + t*s)/2) + this_row = _mat2vec((s*t + t*s)/2) Q_rows.append(W.coordinates(this_row)) Q = matrix(field, W.dimension(), Q_rows) Qs.append(Q) - return Qs + return (Qs, S) def _embed_complex_matrix(M): @@ -986,7 +1179,7 @@ def _embed_complex_matrix(M): """ n = M.nrows() if M.ncols() != n: - raise ArgumentError("the matrix 'M' must be square") + raise ValueError("the matrix 'M' must be square") field = M.base_ring() blocks = [] for z in M.list(): @@ -1014,9 +1207,9 @@ def _unembed_complex_matrix(M): """ n = ZZ(M.nrows()) if M.ncols() != n: - raise ArgumentError("the matrix 'M' must be square") + raise ValueError("the matrix 'M' must be square") if not n.mod(2).is_zero(): - raise ArgumentError("the matrix 'M' must be a complex embedding") + raise ValueError("the matrix 'M' must be a complex embedding") F = QuadraticField(-1, 'i') i = F.gen() @@ -1028,14 +1221,26 @@ def _unembed_complex_matrix(M): for j in xrange(n/2): submat = M[2*k:2*k+2,2*j:2*j+2] if submat[0,0] != submat[1,1]: - raise ArgumentError('bad real submatrix') + raise ValueError('bad real submatrix') if submat[0,1] != -submat[1,0]: - raise ArgumentError('bad imag submatrix') + raise ValueError('bad imag submatrix') z = submat[0,0] + submat[1,0]*i elements.append(z) return matrix(F, n/2, elements) +# The usual inner product on R^n. +def _usual_ip(x,y): + return x.vector().inner_product(y.vector()) + +# The inner product used for the real symmetric simple EJA. +# We keep it as a separate function because e.g. the complex +# algebra uses the same inner product, except divided by 2. +def _matrix_ip(X,Y): + X_mat = X.natural_representation() + Y_mat = Y.natural_representation() + return (X_mat*Y_mat).trace() + def RealSymmetricSimpleEJA(n, field=QQ): """ @@ -1066,9 +1271,13 @@ def RealSymmetricSimpleEJA(n, field=QQ): """ S = _real_symmetric_basis(n, field=field) - Qs = _multiplication_table_from_matrix_basis(S) + (Qs, T) = _multiplication_table_from_matrix_basis(S) - return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=n) + return FiniteDimensionalEuclideanJordanAlgebra(field, + Qs, + rank=n, + natural_basis=T, + inner_product=_matrix_ip) def ComplexHermitianSimpleEJA(n, field=QQ): @@ -1090,8 +1299,22 @@ def ComplexHermitianSimpleEJA(n, field=QQ): """ S = _complex_hermitian_basis(n) - Qs = _multiplication_table_from_matrix_basis(S) - return FiniteDimensionalEuclideanJordanAlgebra(field, Qs, rank=n) + (Qs, T) = _multiplication_table_from_matrix_basis(S) + + # Since a+bi on the diagonal is represented as + # + # a + bi = [ a b ] + # [ -b a ], + # + # we'll double-count the "a" entries if we take the trace of + # the embedding. + ip = lambda X,Y: _matrix_ip(X,Y)/2 + + return FiniteDimensionalEuclideanJordanAlgebra(field, + Qs, + rank=n, + natural_basis=T, + inner_product=ip) def QuaternionHermitianSimpleEJA(n): @@ -1140,10 +1363,10 @@ def JordanSpinSimpleEJA(n, field=QQ): In one dimension, this is the reals under multiplication:: - sage: J1 = JordanSpinSimpleEJA(1) - sage: J2 = eja_rn(1) - sage: J1 == J2 - True + sage: J1 = JordanSpinSimpleEJA(1) + sage: J2 = eja_rn(1) + sage: J1 == J2 + True """ Qs = [] @@ -1162,4 +1385,7 @@ def JordanSpinSimpleEJA(n, field=QQ): # The rank of the spin factor algebra is two, UNLESS we're in a # one-dimensional ambient space (the rank is bounded by the # ambient dimension). - return FiniteDimensionalEuclideanJordanAlgebra(field, Qs, rank=min(n,2)) + return FiniteDimensionalEuclideanJordanAlgebra(field, + Qs, + rank=min(n,2), + inner_product=_usual_ip)