# -*- coding: utf-8 -*- from itertools import izip from sage.matrix.constructor import matrix from sage.modules.free_module import VectorSpace from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement # TODO: make this unnecessary somehow. from sage.misc.lazy_import import lazy_import lazy_import('mjo.eja.eja_algebra', 'FiniteDimensionalEuclideanJordanAlgebra') lazy_import('mjo.eja.eja_element_subalgebra', 'FiniteDimensionalEuclideanJordanElementSubalgebra') from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator from mjo.eja.eja_utils import _mat2vec class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): """ An element of a Euclidean Jordan algebra. """ def __dir__(self): """ Oh man, I should not be doing this. This hides the "disabled" methods ``left_matrix`` and ``matrix`` from introspection; in particular it removes them from tab-completion. """ return filter(lambda s: s not in ['left_matrix', 'matrix'], dir(self.__class__) ) def __pow__(self, n): """ Return ``self`` raised to the power ``n``. Jordan algebras are always power-associative; see for example Faraut and Korányi, Proposition II.1.2 (ii). We have to override this because our superclass uses row vectors instead of column vectors! We, on the other hand, assume column vectors everywhere. SETUP:: sage: from mjo.eja.eja_algebra import random_eja TESTS: The definition of `x^2` is the unambiguous `x*x`:: sage: set_random_seed() sage: x = random_eja().random_element() sage: x*x == (x^2) True A few examples of power-associativity:: sage: set_random_seed() sage: x = random_eja().random_element() sage: x*(x*x)*(x*x) == x^5 True sage: (x*x)*(x*x*x) == x^5 True We also know that powers operator-commute (Koecher, Chapter III, Corollary 1):: sage: set_random_seed() sage: x = random_eja().random_element() sage: m = ZZ.random_element(0,10) sage: n = ZZ.random_element(0,10) sage: Lxm = (x^m).operator() sage: Lxn = (x^n).operator() sage: Lxm*Lxn == Lxn*Lxm True """ if n == 0: return self.parent().one() elif n == 1: return self else: return (self**(n-1))*self def apply_univariate_polynomial(self, p): """ Apply the univariate polynomial ``p`` to this element. A priori, SageMath won't allow us to apply a univariate polynomial to an element of an EJA, because we don't know that EJAs are rings (they are usually not associative). Of course, we know that EJAs are power-associative, so the operation is ultimately kosher. This function sidesteps the CAS to get the answer we want and expect. SETUP:: sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA, ....: random_eja) EXAMPLES:: sage: R = PolynomialRing(QQ, 't') sage: t = R.gen(0) sage: p = t^4 - t^3 + 5*t - 2 sage: J = RealCartesianProductEJA(5) sage: J.one().apply_univariate_polynomial(p) == 3*J.one() True TESTS: We should always get back an element of the algebra:: sage: set_random_seed() sage: p = PolynomialRing(QQ, 't').random_element() sage: J = random_eja() sage: x = J.random_element() sage: x.apply_univariate_polynomial(p) in J True """ if len(p.variables()) > 1: raise ValueError("not a univariate polynomial") P = self.parent() R = P.base_ring() # Convert the coeficcients to the parent's base ring, # because a priori they might live in an (unnecessarily) # larger ring for which P.sum() would fail below. cs = [ R(c) for c in p.coefficients(sparse=False) ] return P.sum( cs[k]*(self**k) for k in range(len(cs)) ) def characteristic_polynomial(self): """ Return the characteristic polynomial of this element. SETUP:: sage: from mjo.eja.eja_algebra import RealCartesianProductEJA EXAMPLES: The rank of `R^3` is three, and the minimal polynomial of the identity element is `(t-1)` from which it follows that the characteristic polynomial should be `(t-1)^3`:: sage: J = RealCartesianProductEJA(3) sage: J.one().characteristic_polynomial() t^3 - 3*t^2 + 3*t - 1 Likewise, the characteristic of the zero element in the rank-three algebra `R^{n}` should be `t^{3}`:: sage: J = RealCartesianProductEJA(3) sage: J.zero().characteristic_polynomial() t^3 TESTS: The characteristic polynomial of an element should evaluate to zero on that element:: sage: set_random_seed() sage: x = RealCartesianProductEJA(3).random_element() sage: p = x.characteristic_polynomial() sage: x.apply_univariate_polynomial(p) 0 The characteristic polynomials of the zero and unit elements should be what we think they are in a subalgebra, too:: sage: J = RealCartesianProductEJA(3) sage: p1 = J.one().characteristic_polynomial() sage: q1 = J.zero().characteristic_polynomial() sage: e0,e1,e2 = J.gens() sage: A = (e0 + 2*e1 + 3*e2).subalgebra_generated_by() # dim 3 sage: p2 = A.one().characteristic_polynomial() sage: q2 = A.zero().characteristic_polynomial() sage: p1 == p2 True sage: q1 == q2 True """ p = self.parent().characteristic_polynomial() return p(*self.to_vector()) def inner_product(self, other): """ Return the parent algebra's inner product of myself and ``other``. SETUP:: sage: from mjo.eja.eja_algebra import ( ....: ComplexHermitianEJA, ....: JordanSpinEJA, ....: QuaternionHermitianEJA, ....: RealSymmetricEJA, ....: random_eja) 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 = JordanSpinEJA(3) sage: x = vector(QQ,[1,2,3]) sage: y = vector(QQ,[4,5,6]) sage: x.inner_product(y) 32 sage: J.from_vector(x).inner_product(J.from_vector(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 = RealSymmetricEJA(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 = ComplexHermitianEJA(3) sage: J.one().inner_product(J.one()) 3 Ditto for the quaternions:: sage: J = QuaternionHermitianEJA(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,y = J.random_elements(2) sage: x.inner_product(y) in RLF 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 with ``other``. SETUP:: sage: from mjo.eja.eja_algebra import random_eja EXAMPLES: The definition of a Jordan algebra says that any element operator-commutes with its square:: sage: set_random_seed() sage: x = random_eja().random_element() sage: x.operator_commutes_with(x^2) True TESTS: Test Lemma 1 from Chapter III of Koecher:: sage: set_random_seed() sage: u,v = random_eja().random_elements(2) sage: lhs = u.operator_commutes_with(u*v) sage: rhs = v.operator_commutes_with(u^2) sage: lhs == rhs True Test the first polarization identity from my notes, Koecher Chapter III, or from Baes (2.3):: sage: set_random_seed() sage: x,y = random_eja().random_elements(2) sage: Lx = x.operator() sage: Ly = y.operator() sage: Lxx = (x*x).operator() sage: Lxy = (x*y).operator() sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly) True Test the second polarization identity from my notes or from Baes (2.4):: sage: set_random_seed() sage: x,y,z = random_eja().random_elements(3) sage: Lx = x.operator() sage: Ly = y.operator() sage: Lz = z.operator() sage: Lzy = (z*y).operator() sage: Lxy = (x*y).operator() sage: Lxz = (x*z).operator() sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly) True Test the third polarization identity from my notes or from Baes (2.5):: sage: set_random_seed() sage: u,y,z = random_eja().random_elements(3) sage: Lu = u.operator() sage: Ly = y.operator() sage: Lz = z.operator() sage: Lzy = (z*y).operator() sage: Luy = (u*y).operator() sage: Luz = (u*z).operator() sage: Luyz = (u*(y*z)).operator() sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly sage: bool(lhs == rhs) True """ if not other in self.parent(): raise TypeError("'other' must live in the same algebra") A = self.operator() B = other.operator() return (A*B == B*A) def det(self): """ Return my determinant, the product of my eigenvalues. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: random_eja) EXAMPLES:: sage: J = JordanSpinEJA(2) sage: e0,e1 = J.gens() sage: x = sum( J.gens() ) sage: x.det() 0 :: sage: J = JordanSpinEJA(3) sage: e0,e1,e2 = J.gens() sage: x = sum( J.gens() ) sage: x.det() -1 TESTS: An element is invertible if and only if its determinant is non-zero:: sage: set_random_seed() sage: x = random_eja().random_element() sage: x.is_invertible() == (x.det() != 0) True Ensure that the determinant is multiplicative on an associative subalgebra as in Faraut and Korányi's Proposition II.2.2:: sage: set_random_seed() sage: J = random_eja().random_element().subalgebra_generated_by() sage: x,y = J.random_elements(2) sage: (x*y).det() == x.det()*y.det() True """ P = self.parent() r = P.rank() p = P._charpoly_coeff(0) # The _charpoly_coeff function already adds the factor of # -1 to ensure that _charpoly_coeff(0) is really what # appears in front of t^{0} in the charpoly. However, # we want (-1)^r times THAT for the determinant. return ((-1)**r)*p(*self.to_vector()) def inverse(self): """ Return the Jordan-multiplicative inverse of this element. ALGORITHM: We appeal to the quadratic representation as in Koecher's Theorem 12 in Chapter III, Section 5. SETUP:: sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, ....: JordanSpinEJA, ....: random_eja) EXAMPLES: The inverse in the spin factor algebra is given in Alizadeh's Example 11.11:: sage: set_random_seed() sage: J = JordanSpinEJA.random_instance() sage: x = J.random_element() sage: while not x.is_invertible(): ....: x = J.random_element() sage: x_vec = x.to_vector() sage: x0 = x_vec[0] sage: x_bar = x_vec[1:] sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar)) sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list()) sage: x_inverse = coeff*inv_vec sage: x.inverse() == J.from_vector(x_inverse) True Trying to invert a non-invertible element throws an error: sage: JordanSpinEJA(3).zero().inverse() Traceback (most recent call last): ... ValueError: element is not invertible TESTS: The identity element is its own inverse:: sage: set_random_seed() sage: J = random_eja() sage: J.one().inverse() == J.one() True If an element has an inverse, it acts like one:: sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: (not x.is_invertible()) or (x.inverse()*x == J.one()) True The inverse of the inverse is what we started with:: sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: (not x.is_invertible()) or (x.inverse().inverse() == x) True Proposition II.2.3 in Faraut and Korányi says that the inverse of an element is the inverse of its left-multiplication operator applied to the algebra's identity, when that inverse exists:: sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: (not x.operator().is_invertible()) or ( ....: x.operator().inverse()(J.one()) == x.inverse() ) True Proposition II.2.4 in Faraut and Korányi gives a formula for the inverse based on the characteristic polynomial and the Cayley-Hamilton theorem for Euclidean Jordan algebras:: sage: set_random_seed() sage: J = ComplexHermitianEJA(3) sage: x = J.random_element() sage: while not x.is_invertible(): ....: x = J.random_element() sage: r = J.rank() sage: a = x.characteristic_polynomial().coefficients(sparse=False) sage: expected = (-1)^(r+1)/x.det() sage: expected *= sum( a[i+1]*x^i for i in range(r) ) sage: x.inverse() == expected True """ if not self.is_invertible(): raise ValueError("element is not invertible") return (~self.quadratic_representation())(self) def is_invertible(self): """ Return whether or not this element is invertible. ALGORITHM: The usual way to do this is to check if the determinant is zero, but we need the characteristic polynomial for the determinant. The minimal polynomial is a lot easier to get, so we use Corollary 2 in Chapter V of Koecher to check whether or not the paren't algebra's zero element is a root of this element's minimal polynomial. Beware that we can't use the superclass method, because it relies on the algebra being associative. SETUP:: sage: from mjo.eja.eja_algebra import random_eja TESTS: The identity element is always invertible:: sage: set_random_seed() sage: J = random_eja() sage: J.one().is_invertible() True The zero element is never invertible in a non-trivial algebra:: sage: set_random_seed() sage: J = random_eja() sage: (not J.is_trivial()) and J.zero().is_invertible() False """ if self.is_zero(): if self.parent().is_trivial(): return True else: return False # In fact, we only need to know if the constant term is non-zero, # so we can pass in the field's zero element instead. zero = self.base_ring().zero() p = self.minimal_polynomial() return not (p(zero) == zero) def is_primitive_idempotent(self): """ Return whether or not this element is a primitive (or minimal) idempotent. A primitive idempotent is a non-zero idempotent that is not the sum of two other non-zero idempotents. Remark 2.7.15 in Baes shows that this is what he refers to as a "minimal idempotent." An element of a Euclidean Jordan algebra is a minimal idempotent if it :meth:`is_idempotent` and if its Peirce subalgebra corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes, Proposition 2.7.17). SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: RealSymmetricEJA, ....: TrivialEJA, ....: random_eja) WARNING:: This method is sloooooow. EXAMPLES: The spectral decomposition of a non-regular element should always contain at least one non-minimal idempotent:: sage: J = RealSymmetricEJA(3, AA) sage: x = sum(J.gens()) sage: x.is_regular() False sage: [ c.is_primitive_idempotent() ....: for (l,c) in x.spectral_decomposition() ] [False, True] On the other hand, the spectral decomposition of a regular element should always be in terms of minimal idempotents:: sage: J = JordanSpinEJA(4, AA) sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) ) sage: x.is_regular() True sage: [ c.is_primitive_idempotent() ....: for (l,c) in x.spectral_decomposition() ] [True, True] TESTS: The identity element is minimal only in an EJA of rank one:: sage: set_random_seed() sage: J = random_eja() sage: J.rank() == 1 or not J.one().is_primitive_idempotent() True A non-idempotent cannot be a minimal idempotent:: sage: set_random_seed() sage: J = JordanSpinEJA(4) sage: x = J.random_element() sage: (not x.is_idempotent()) and x.is_primitive_idempotent() False Proposition 2.7.19 in Baes says that an element is a minimal idempotent if and only if it's idempotent with trace equal to unity:: sage: set_random_seed() sage: J = JordanSpinEJA(4) sage: x = J.random_element() sage: expected = (x.is_idempotent() and x.trace() == 1) sage: actual = x.is_primitive_idempotent() sage: actual == expected True Primitive idempotents must be non-zero:: sage: set_random_seed() sage: J = random_eja() sage: J.zero().is_idempotent() True sage: J.zero().is_primitive_idempotent() False As a consequence of the fact that primitive idempotents must be non-zero, there are no primitive idempotents in a trivial Euclidean Jordan algebra:: sage: J = TrivialEJA() sage: J.one().is_idempotent() True sage: J.one().is_primitive_idempotent() False """ if not self.is_idempotent(): return False if self.is_zero(): return False (_,_,J1) = self.parent().peirce_decomposition(self) return (J1.dimension() == 1) def is_nilpotent(self): """ Return whether or not some power of this element is zero. ALGORITHM: We use Theorem 5 in Chapter III of Koecher, which says that an element ``x`` is nilpotent if and only if ``x.operator()`` is nilpotent. And it is a basic fact of linear algebra that an operator on an `n`-dimensional space is nilpotent if and only if, when raised to the `n`th power, it equals the zero operator (for example, see Axler Corollary 8.8). SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: random_eja) EXAMPLES:: sage: J = JordanSpinEJA(3) sage: x = sum(J.gens()) sage: x.is_nilpotent() False TESTS: The identity element is never nilpotent, except in a trivial EJA:: sage: set_random_seed() sage: J = random_eja() sage: J.one().is_nilpotent() and not J.is_trivial() False The additive identity is always nilpotent:: sage: set_random_seed() sage: random_eja().zero().is_nilpotent() True """ P = self.parent() zero_operator = P.zero().operator() return self.operator()**P.dimension() == zero_operator def is_regular(self): """ Return whether or not this is a regular element. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: random_eja) EXAMPLES: The identity element always has degree one, but any element linearly-independent from it is regular:: sage: J = JordanSpinEJA(5) sage: J.one().is_regular() False sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity sage: for x in J.gens(): ....: (J.one() + x).is_regular() False True True True True TESTS: The zero element should never be regular, unless the parent algebra has dimension less than or equal to one:: sage: set_random_seed() sage: J = random_eja() sage: J.dimension() <= 1 or not J.zero().is_regular() True The unit element isn't regular unless the algebra happens to consist of only its scalar multiples:: sage: set_random_seed() sage: J = random_eja() sage: J.dimension() <= 1 or not J.one().is_regular() True """ return self.degree() == self.parent().rank() def degree(self): """ Return the degree of this element, which is defined to be the degree of its minimal polynomial. ALGORITHM: For now, we skip the messy minimal polynomial computation and instead return the dimension of the vector space spanned by the powers of this element. The latter is a bit more straightforward to compute. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: random_eja) EXAMPLES:: sage: J = JordanSpinEJA(4) sage: J.one().degree() 1 sage: e0,e1,e2,e3 = J.gens() sage: (e0 - e1).degree() 2 In the spin factor algebra (of rank two), all elements that aren't multiples of the identity are regular:: sage: set_random_seed() sage: J = JordanSpinEJA.random_instance() sage: x = J.random_element() sage: x == x.coefficient(0)*J.one() or x.degree() == 2 True TESTS: The zero and unit elements are both of degree one in nontrivial algebras:: sage: set_random_seed() sage: J = random_eja() sage: d = J.zero().degree() sage: (J.is_trivial() and d == 0) or d == 1 True sage: d = J.one().degree() sage: (J.is_trivial() and d == 0) or d == 1 True Our implementation agrees with the definition:: sage: set_random_seed() sage: x = random_eja().random_element() sage: x.degree() == x.minimal_polynomial().degree() True """ if self.is_zero() and not self.parent().is_trivial(): # The minimal polynomial of zero in a nontrivial algebra # is "t"; in a trivial algebra it's "1" by convention # (it's an empty product). return 1 return self.subalgebra_generated_by().dimension() def left_matrix(self): """ Our parent class defines ``left_matrix`` and ``matrix`` methods whose names are misleading. We don't want them. """ raise NotImplementedError("use operator().matrix() instead") matrix = left_matrix def minimal_polynomial(self): """ Return the minimal polynomial of this element, as a function of the variable `t`. ALGORITHM: We restrict ourselves to the associative subalgebra generated by this element, and then return the minimal polynomial of this element's operator matrix (in that subalgebra). This works by Baes Proposition 2.3.16. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: RealSymmetricEJA, ....: TrivialEJA, ....: random_eja) EXAMPLES: Keeping in mind that the polynomial ``1`` evaluates the identity element (also the zero element) of the trivial algebra, it is clear that the polynomial ``1`` is the minimal polynomial of the only element in a trivial algebra:: sage: J = TrivialEJA() sage: J.one().minimal_polynomial() 1 sage: J.zero().minimal_polynomial() 1 TESTS: The minimal polynomial of the identity and zero elements are always the same:: sage: set_random_seed() sage: J = random_eja(nontrivial=True) sage: J.one().minimal_polynomial() t - 1 sage: J.zero().minimal_polynomial() t The degree of an element is (by one definition) the degree of its minimal polynomial:: 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. We require the dimension of the algebra to be at least two here so that said elements actually exist:: sage: set_random_seed() sage: n_max = max(2, JordanSpinEJA._max_test_case_size()) sage: n = ZZ.random_element(2, n_max) sage: J = JordanSpinEJA(n) sage: y = J.random_element() sage: while y == y.coefficient(0)*J.one(): ....: y = J.random_element() sage: y0 = y.to_vector()[0] sage: y_bar = y.to_vector()[1:] sage: actual = y.minimal_polynomial() sage: t = PolynomialRing(J.base_ring(),'t').gen(0) sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2) sage: bool(actual == expected) True The minimal polynomial should always kill its element:: sage: set_random_seed() sage: x = random_eja().random_element() sage: p = x.minimal_polynomial() sage: x.apply_univariate_polynomial(p) 0 The minimal polynomial is invariant under a change of basis, and in particular, a re-scaling of the basis:: sage: set_random_seed() sage: n_max = RealSymmetricEJA._max_test_case_size() sage: n = ZZ.random_element(1, n_max) sage: J1 = RealSymmetricEJA(n,QQ) sage: J2 = RealSymmetricEJA(n,QQ,normalize_basis=False) sage: X = random_matrix(QQ,n) sage: X = X*X.transpose() sage: x1 = J1(X) sage: x2 = J2(X) sage: x1.minimal_polynomial() == x2.minimal_polynomial() True """ if self.is_zero(): # We would generate a zero-dimensional subalgebra # where the minimal polynomial would be constant. # That might be correct, but only if *this* algebra # is trivial too. if not self.parent().is_trivial(): # Pretty sure we know what the minimal polynomial of # the zero operator is going to be. This ensures # consistency of e.g. the polynomial variable returned # in the "normal" case without us having to think about it. return self.operator().minimal_polynomial() A = self.subalgebra_generated_by() return A(self).operator().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. SETUP:: sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, ....: QuaternionHermitianEJA) EXAMPLES:: sage: J = ComplexHermitianEJA(3) sage: J.one() e0 + e3 + 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] :: sage: J = QuaternionHermitianEJA(3) sage: J.one() e0 + e5 + e14 sage: J.one().natural_representation() [1 0 0 0 0 0 0 0 0 0 0 0] [0 1 0 0 0 0 0 0 0 0 0 0] [0 0 1 0 0 0 0 0 0 0 0 0] [0 0 0 1 0 0 0 0 0 0 0 0] [0 0 0 0 1 0 0 0 0 0 0 0] [0 0 0 0 0 1 0 0 0 0 0 0] [0 0 0 0 0 0 1 0 0 0 0 0] [0 0 0 0 0 0 0 1 0 0 0 0] [0 0 0 0 0 0 0 0 1 0 0 0] [0 0 0 0 0 0 0 0 0 1 0 0] [0 0 0 0 0 0 0 0 0 0 1 0] [0 0 0 0 0 0 0 0 0 0 0 1] """ B = self.parent().natural_basis() W = self.parent().natural_basis_space() return W.linear_combination(izip(B,self.to_vector())) def norm(self): """ The norm of this element with respect to :meth:`inner_product`. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: RealCartesianProductEJA) EXAMPLES:: sage: J = RealCartesianProductEJA(2) sage: x = sum(J.gens()) sage: x.norm() sqrt(2) :: sage: J = JordanSpinEJA(4) sage: x = sum(J.gens()) sage: x.norm() 2 """ return self.inner_product(self).sqrt() def operator(self): """ Return the left-multiplication-by-this-element operator on the ambient algebra. SETUP:: sage: from mjo.eja.eja_algebra import random_eja TESTS:: sage: set_random_seed() sage: J = random_eja() sage: x,y = J.random_elements(2) sage: x.operator()(y) == x*y True sage: y.operator()(x) == x*y True """ P = self.parent() left_mult_by_self = lambda y: self*y L = P.module_morphism(function=left_mult_by_self, codomain=P) return FiniteDimensionalEuclideanJordanAlgebraOperator( P, P, L.matrix() ) def quadratic_representation(self, other=None): """ Return the quadratic representation of this element. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: random_eja) EXAMPLES: The explicit form in the spin factor algebra is given by Alizadeh's Example 11.12:: sage: set_random_seed() sage: x = JordanSpinEJA.random_instance().random_element() sage: x_vec = x.to_vector() sage: n = x_vec.degree() sage: x0 = x_vec[0] sage: x_bar = x_vec[1:] sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)]) sage: B = 2*x0*x_bar.row() sage: C = 2*x0*x_bar.column() sage: D = matrix.identity(QQ, n-1) sage: D = (x0^2 - x_bar.inner_product(x_bar))*D sage: D = D + 2*x_bar.tensor_product(x_bar) sage: Q = matrix.block(2,2,[A,B,C,D]) sage: Q == x.quadratic_representation().matrix() True Test all of the properties from Theorem 11.2 in Alizadeh:: sage: set_random_seed() sage: J = random_eja() sage: x,y = J.random_elements(2) sage: Lx = x.operator() sage: Lxx = (x*x).operator() sage: Qx = x.quadratic_representation() sage: Qy = y.quadratic_representation() sage: Qxy = x.quadratic_representation(y) sage: Qex = J.one().quadratic_representation(x) sage: n = ZZ.random_element(10) sage: Qxn = (x^n).quadratic_representation() Property 1: sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy True Property 2 (multiply on the right for :trac:`28272`): sage: alpha = J.base_ring().random_element() sage: (alpha*x).quadratic_representation() == Qx*(alpha^2) True Property 3: sage: not x.is_invertible() or ( Qx(x.inverse()) == x ) True sage: not x.is_invertible() or ( ....: ~Qx ....: == ....: x.inverse().quadratic_representation() ) True sage: Qxy(J.one()) == x*y True Property 4: sage: not x.is_invertible() or ( ....: x.quadratic_representation(x.inverse())*Qx ....: == Qx*x.quadratic_representation(x.inverse()) ) True sage: not x.is_invertible() or ( ....: x.quadratic_representation(x.inverse())*Qx ....: == ....: 2*Lx*Qex - Qx ) True sage: 2*Lx*Qex - Qx == Lxx True Property 5: sage: Qy(x).quadratic_representation() == Qy*Qx*Qy True Property 6: sage: Qxn == (Qx)^n True Property 7: sage: not x.is_invertible() or ( ....: Qx*x.inverse().operator() == Lx ) True Property 8: sage: not x.operator_commutes_with(y) or ( ....: Qx(y)^n == Qxn(y^n) ) True """ if other is None: other=self elif not other in self.parent(): raise TypeError("'other' must live in the same algebra") L = self.operator() M = other.operator() return ( L*M + M*L - (self*other).operator() ) def spectral_decomposition(self): """ Return the unique spectral decomposition of this element. ALGORITHM: Following Faraut and Korányi's Theorem III.1.1, we restrict this element's left-multiplication-by operator to the subalgebra it generates. We then compute the spectral decomposition of that operator, and the spectral projectors we get back must be the left-multiplication-by operators for the idempotents we seek. Thus applying them to the identity element gives us those idempotents. Since the eigenvalues are required to be distinct, we take the spectral decomposition of the zero element to be zero times the identity element of the algebra (which is idempotent, obviously). SETUP:: sage: from mjo.eja.eja_algebra import RealSymmetricEJA EXAMPLES: The spectral decomposition of the identity is ``1`` times itself, and the spectral decomposition of zero is ``0`` times the identity:: sage: J = RealSymmetricEJA(3,AA) sage: J.one() e0 + e2 + e5 sage: J.one().spectral_decomposition() [(1, e0 + e2 + e5)] sage: J.zero().spectral_decomposition() [(0, e0 + e2 + e5)] TESTS:: sage: J = RealSymmetricEJA(4,AA) sage: x = sum(J.gens()) sage: sd = x.spectral_decomposition() sage: l0 = sd[0][0] sage: l1 = sd[1][0] sage: c0 = sd[0][1] sage: c1 = sd[1][1] sage: c0.inner_product(c1) == 0 True sage: c0.is_idempotent() True sage: c1.is_idempotent() True sage: c0 + c1 == J.one() True sage: l0*c0 + l1*c1 == x True """ P = self.parent() A = self.subalgebra_generated_by(orthonormalize_basis=True) result = [] for (evalue, proj) in A(self).operator().spectral_decomposition(): result.append( (evalue, proj(A.one()).superalgebra_element()) ) return result def subalgebra_generated_by(self, orthonormalize_basis=False): """ Return the associative subalgebra of the parent EJA generated by this element. Since our parent algebra is unital, we want "subalgebra" to mean "unital subalgebra" as well; thus the subalgebra that an element generates will itself be a Euclidean Jordan algebra after restricting the algebra operations appropriately. This is the subalgebra that Faraut and Korányi work with in section II.2, for example. SETUP:: sage: from mjo.eja.eja_algebra import random_eja TESTS: This subalgebra, being composed of only powers, is associative:: sage: set_random_seed() sage: x0 = random_eja().random_element() sage: A = x0.subalgebra_generated_by() sage: x,y,z = A.random_elements(3) sage: (x*y)*z == x*(y*z) True Squaring in the subalgebra should work the same as in the superalgebra:: sage: set_random_seed() sage: x = random_eja().random_element() sage: A = x.subalgebra_generated_by() sage: A(x^2) == A(x)*A(x) True By definition, the subalgebra generated by the zero element is the one-dimensional algebra generated by the identity element... unless the original algebra was trivial, in which case the subalgebra is trivial too:: sage: set_random_seed() sage: A = random_eja().zero().subalgebra_generated_by() sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1 True """ return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis) def subalgebra_idempotent(self): """ Find an idempotent in the associative subalgebra I generate using Proposition 2.3.5 in Baes. SETUP:: sage: from mjo.eja.eja_algebra import random_eja TESTS: Ensure that we can find an idempotent in a non-trivial algebra where there are non-nilpotent elements:: sage: set_random_seed() sage: J = random_eja(nontrivial=True) sage: x = J.random_element() sage: while x.is_nilpotent(): ....: x = J.random_element() sage: c = x.subalgebra_idempotent() sage: c^2 == c True """ if self.parent().is_trivial(): return self if self.is_nilpotent(): raise ValueError("this only works with non-nilpotent elements!") J = self.subalgebra_generated_by() u = J(self) # The image of the matrix of left-u^m-multiplication # will be minimal for some natural number s... s = 0 minimal_dim = J.dimension() for i in xrange(1, minimal_dim): this_dim = (u**i).operator().matrix().image().dimension() if this_dim < minimal_dim: minimal_dim = this_dim s = i # Now minimal_matrix should correspond to the smallest # non-zero subspace in Baes's (or really, Koecher's) # proposition. # # However, we need to restrict the matrix to work on the # subspace... or do we? Can't we just solve, knowing that # A(c) = u^(s+1) should have a solution in the big space, # too? # # Beware, solve_right() means that we're using COLUMN vectors. # Our FiniteDimensionalAlgebraElement superclass uses rows. u_next = u**(s+1) A = u_next.operator().matrix() c = J.from_vector(A.solve_right(u_next.to_vector())) # Now c is the idempotent we want, but it still lives in the subalgebra. return c.superalgebra_element() def trace(self): """ Return my trace, the sum of my eigenvalues. In a trivial algebra, however you want to look at it, the trace is an empty sum for which we declare the result to be zero. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: RealCartesianProductEJA, ....: TrivialEJA, ....: random_eja) EXAMPLES:: sage: J = TrivialEJA() sage: J.zero().trace() 0 :: sage: J = JordanSpinEJA(3) sage: x = sum(J.gens()) sage: x.trace() 2 :: sage: J = RealCartesianProductEJA(5) sage: J.one().trace() 5 TESTS: The trace of an element is a real number:: sage: set_random_seed() sage: J = random_eja() sage: J.random_element().trace() in RLF True """ P = self.parent() r = P.rank() if r == 0: # Special case for the trivial algebra where # the trace is an empty sum. return P.base_ring().zero() p = P._charpoly_coeff(r-1) # The _charpoly_coeff function already adds the factor of # -1 to ensure that _charpoly_coeff(r-1) is really what # appears in front of t^{r-1} in the charpoly. However, # we want the negative of THAT for the trace. return -p(*self.to_vector()) def trace_inner_product(self, other): """ Return the trace inner product of myself and ``other``. SETUP:: sage: from mjo.eja.eja_algebra import random_eja TESTS: The trace inner product is commutative, bilinear, and associative:: sage: set_random_seed() sage: J = random_eja() sage: x,y,z = J.random_elements(3) sage: # commutative sage: x.trace_inner_product(y) == y.trace_inner_product(x) True sage: # bilinear sage: a = J.base_ring().random_element(); sage: actual = (a*(x+z)).trace_inner_product(y) sage: expected = ( a*x.trace_inner_product(y) + ....: a*z.trace_inner_product(y) ) sage: actual == expected True sage: actual = x.trace_inner_product(a*(y+z)) sage: expected = ( a*x.trace_inner_product(y) + ....: a*x.trace_inner_product(z) ) sage: actual == expected True sage: # associative sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z) True """ if not other in self.parent(): raise TypeError("'other' must live in the same algebra") return (self*other).trace() def trace_norm(self): """ The norm of this element with respect to :meth:`trace_inner_product`. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: RealCartesianProductEJA) EXAMPLES:: sage: J = RealCartesianProductEJA(2) sage: x = sum(J.gens()) sage: x.trace_norm() sqrt(2) :: sage: J = JordanSpinEJA(4) sage: x = sum(J.gens()) sage: x.trace_norm() 2*sqrt(2) """ return self.trace_inner_product(self).sqrt()