From f9e1e977040c2c3b1019a0bcc5776384a13e96c4 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 29 Jul 2019 20:31:33 -0400 Subject: [PATCH] eja: factor out the EJA element class into its own module. --- mjo/eja/eja_algebra.py | 1249 +--------------------------------------- mjo/eja/eja_element.py | 1240 +++++++++++++++++++++++++++++++++++++++ mjo/eja/eja_utils.py | 9 + 3 files changed, 1254 insertions(+), 1244 deletions(-) create mode 100644 mjo/eja/eja_element.py create mode 100644 mjo/eja/eja_utils.py diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index ad2afbd..fb840ed 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -5,18 +5,13 @@ are used in optimization, and have some additional nice methods beyond what can be supported in a general Jordan Algebra. """ - - from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra -from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element import FiniteDimensionalAlgebraElement from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra from sage.categories.finite_dimensional_algebras_with_basis import FiniteDimensionalAlgebrasWithBasis -from sage.functions.other import sqrt from sage.matrix.constructor import matrix from sage.misc.cachefunc import cached_method from sage.misc.prandom import choice from sage.modules.free_module import VectorSpace -from sage.modules.free_module_element import vector from sage.rings.integer_ring import ZZ from sage.rings.number_field.number_field import QuadraticField from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing @@ -24,8 +19,8 @@ from sage.rings.rational_field import QQ from sage.structure.element import is_Matrix from sage.structure.category_object import normalize_names -from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator - +from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement +from mjo.eja.eja_utils import _vec2mat, _mat2vec class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): @staticmethod @@ -56,7 +51,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): return fda.__classcall__(cls, field, mult_table, - rank=rank, + rank, assume_associative=assume_associative, names=names, category=cat, @@ -214,7 +209,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): R = PolynomialRing(self.base_ring(), names) J = FiniteDimensionalEuclideanJordanAlgebra(R, self._multiplication_table, - rank=r) + r) idmat = matrix.identity(J.base_ring(), n) @@ -461,1236 +456,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): return self.zero().vector().parent().ambient_vector_space() - class Element(FiniteDimensionalAlgebraElement): - """ - 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 __init__(self, A, elt=None): - """ - - SETUP:: - - sage: from mjo.eja.eja_algebra import (RealSymmetricEJA, - ....: random_eja) - - EXAMPLES: - - The identity in `S^n` is converted to the identity in the EJA:: - - sage: J = RealSymmetricEJA(3) - sage: I = matrix.identity(QQ,3) - sage: J(I) == J.one() - True - - This skew-symmetric matrix can't be represented in the EJA:: - - sage: J = RealSymmetricEJA(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 - - TESTS: - - Ensure that we can convert any element of the parent's - underlying vector space back into an algebra element whose - vector representation is what we started with:: - - sage: set_random_seed() - sage: J = random_eja() - sage: v = J.vector_space().random_element() - sage: J(v).vector() == v - True - - """ - # 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``. - - Jordan algebras are always power-associative; see for - example Faraut and Koranyi, 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.operator()**(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 - - """ - p = self.parent().characteristic_polynomial() - return p(*self.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(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 = 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 = 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 - 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: J = random_eja() - sage: u = J.random_element() - sage: v = J.random_element() - 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: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - 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: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - 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: J = random_eja() - sage: u = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - 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 - - """ - 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.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 (JordanSpinEJA, - ....: random_eja) - - EXAMPLES: - - The inverse in the spin factor algebra is given in Alizadeh's - Example 11.11:: - - sage: set_random_seed() - sage: n = ZZ.random_element(1,10) - sage: J = JordanSpinEJA(n) - sage: x = J.random_element() - sage: while not x.is_invertible(): - ....: x = J.random_element() - sage: x_vec = x.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(x_inverse) - True - - 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 - - The zero element is never invertible:: - - sage: set_random_seed() - sage: J = random_eja().zero().inverse() - Traceback (most recent call last): - ... - ValueError: element is not invertible - - """ - 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:: - - sage: set_random_seed() - sage: J = random_eja() - sage: J.zero().is_invertible() - False - - """ - zero = self.parent().zero() - p = self.minimal_polynomial() - return not (p(zero) == zero) - - - 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:: - - sage: set_random_seed() - sage: random_eja().one().is_nilpotent() - 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:: - - sage: set_random_seed() - sage: J = random_eja() - sage: J.zero().is_regular() - False - - 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: n = ZZ.random_element(1,10) - sage: J = JordanSpinEJA(n) - 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:: - - sage: set_random_seed() - sage: J = random_eja() - sage: J.zero().degree() - 1 - sage: J.one().degree() - 1 - - Our implementation agrees with the definition:: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: x.degree() == x.minimal_polynomial().degree() - True - - """ - return self.span_of_powers().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, - ....: random_eja) - - TESTS: - - The minimal polynomial of the identity and zero elements are - always the same:: - - sage: set_random_seed() - sage: J = random_eja() - 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:: - - sage: set_random_seed() - sage: n = ZZ.random_element(2,10) - sage: J = JordanSpinEJA(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: 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 - - """ - 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())) - return elt.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 + 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] - - :: - - sage: J = QuaternionHermitianEJA(3) - sage: J.one() - e0 + e9 + 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 = B[0].matrix_space() - return W.linear_combination(zip(self.vector(), B)) - - - 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 = J.random_element() - sage: y = J.random_element() - sage: x.operator()(y) == x*y - True - sage: y.operator()(x) == x*y - True - - """ - P = self.parent() - fda_elt = FiniteDimensionalAlgebraElement(P, self) - return FiniteDimensionalEuclideanJordanAlgebraOperator( - P, - P, - fda_elt.matrix().transpose() ) - - - 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: n = ZZ.random_element(1,10) - sage: J = JordanSpinEJA(n) - sage: x = J.random_element() - sage: x_vec = x.vector() - 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 = J.random_element() - sage: y = J.random_element() - 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 = QQ.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*x.operator()*Qex - Qx ) - True - - sage: 2*x.operator()*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 span_of_powers(self): - """ - Return the vector space spanned by successive powers of - this element. - """ - # The dimension of the subalgebra can't be greater than - # the big algebra, so just put everything into a list - # and let span() get rid of the excess. - # - # We do the extra ambient_vector_space() in case we're messing - # with polynomials and the direct parent is a module. - V = self.parent().vector_space() - return V.span( (self**d).vector() for d in xrange(V.dimension()) ) - - - def subalgebra_generated_by(self): - """ - Return the associative subalgebra of the parent EJA generated - by this element. - - SETUP:: - - sage: from mjo.eja.eja_algebra import random_eja - - TESTS:: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: x.subalgebra_generated_by().is_associative() - True - - Squaring in the subalgebra should work the same as in - the superalgebra:: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: u = x.subalgebra_generated_by().random_element() - sage: u.operator()(u) == u^2 - True - - """ - # First get the subspace spanned by the powers of myself... - V = self.span_of_powers() - F = self.base_ring() - - # Now figure out the entries of the right-multiplication - # matrix for the successive basis elements b0, b1,... of - # that subspace. - mats = [] - for b_right in V.basis(): - eja_b_right = self.parent()(b_right) - b_right_rows = [] - # The first row of the right-multiplication matrix by - # b1 is what we get if we apply that matrix to b1. The - # second row of the right multiplication matrix by b1 - # is what we get when we apply that matrix to b2... - # - # IMPORTANT: this assumes that all vectors are COLUMN - # vectors, unlike our superclass (which uses row vectors). - for b_left in V.basis(): - eja_b_left = self.parent()(b_left) - # Multiply in the original EJA, but then get the - # coordinates from the subalgebra in terms of its - # basis. - this_row = V.coordinates((eja_b_left*eja_b_right).vector()) - b_right_rows.append(this_row) - b_right_matrix = matrix(F, b_right_rows) - mats.append(b_right_matrix) - - # It's an algebra of polynomials in one element, and EJAs - # are power-associative. - # - # TODO: choose generator names intelligently. - # - # The rank is the highest possible degree of a minimal polynomial, - # and is bounded above by the dimension. We know in this case that - # there's an element whose minimal polynomial has the same degree - # as the space's dimension, so that must be its rank too. - return FiniteDimensionalEuclideanJordanAlgebra( - F, - mats, - V.dimension(), - assume_associative=True, - names='f') - - - 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:: - - sage: set_random_seed() - sage: J = random_eja() - 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.is_nilpotent(): - raise ValueError("this only works with non-nilpotent elements!") - - V = self.span_of_powers() - J = 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! - u = J(V.coordinates(self.vector())) - - # The image of the matrix of left-u^m-multiplication - # will be minimal for some natural number s... - s = 0 - minimal_dim = V.dimension() - for i in xrange(1, V.dimension()): - 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_coordinates = A.solve_right(u_next.vector()) - - # Now c_coordinates is the idempotent we want, but it's in - # the coordinate system of the subalgebra. - # - # We need the basis for J, but as elements of the parent algebra. - # - basis = [self.parent(v) for v in V.basis()] - return self.parent().linear_combination(zip(c_coordinates, basis)) - - - def trace(self): - """ - Return my trace, the sum of my eigenvalues. - - SETUP:: - - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, - ....: RealCartesianProductEJA, - ....: random_eja) - - EXAMPLES:: - - 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 J.base_ring() - True - - """ - P = self.parent() - r = P.rank() - 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.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:: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element(); y = J.random_element() - sage: x.trace_inner_product(y) == y.trace_inner_product(x) - True - - The trace inner product is bilinear:: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - sage: a = QQ.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 - - The trace inner product satisfies the compatibility - condition in the definition of a Euclidean Jordan algebra:: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - 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() + Element = FiniteDimensionalEuclideanJordanAlgebraElement class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra): @@ -1900,11 +666,6 @@ def _quaternion_hermitian_basis(n, field=QQ): 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): """ diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py new file mode 100644 index 0000000..1b5e314 --- /dev/null +++ b/mjo/eja/eja_element.py @@ -0,0 +1,1240 @@ +from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element import FiniteDimensionalAlgebraElement +from sage.matrix.constructor import matrix +from sage.modules.free_module import VectorSpace + +# TODO: make this unnecessary somehow. +from sage.misc.lazy_import import lazy_import +lazy_import('mjo.eja.eja_algebra', 'FiniteDimensionalEuclideanJordanAlgebra') +from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator +from mjo.eja.eja_utils import _mat2vec + +class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraElement): + """ + 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 __init__(self, A, elt=None): + """ + + SETUP:: + + sage: from mjo.eja.eja_algebra import (RealSymmetricEJA, + ....: random_eja) + + EXAMPLES: + + The identity in `S^n` is converted to the identity in the EJA:: + + sage: J = RealSymmetricEJA(3) + sage: I = matrix.identity(QQ,3) + sage: J(I) == J.one() + True + + This skew-symmetric matrix can't be represented in the EJA:: + + sage: J = RealSymmetricEJA(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 + + TESTS: + + Ensure that we can convert any element of the parent's + underlying vector space back into an algebra element whose + vector representation is what we started with:: + + sage: set_random_seed() + sage: J = random_eja() + sage: v = J.vector_space().random_element() + sage: J(v).vector() == v + True + + """ + # 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``. + + Jordan algebras are always power-associative; see for + example Faraut and Koranyi, 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.operator()**(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 + + """ + p = self.parent().characteristic_polynomial() + return p(*self.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(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 = 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 = 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 + 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: J = random_eja() + sage: u = J.random_element() + sage: v = J.random_element() + 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: J = random_eja() + sage: x = J.random_element() + sage: y = J.random_element() + 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: J = random_eja() + sage: x = J.random_element() + sage: y = J.random_element() + sage: z = J.random_element() + 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: J = random_eja() + sage: u = J.random_element() + sage: y = J.random_element() + sage: z = J.random_element() + 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 + + """ + 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.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 (JordanSpinEJA, + ....: random_eja) + + EXAMPLES: + + The inverse in the spin factor algebra is given in Alizadeh's + Example 11.11:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,10) + sage: J = JordanSpinEJA(n) + sage: x = J.random_element() + sage: while not x.is_invertible(): + ....: x = J.random_element() + sage: x_vec = x.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(x_inverse) + True + + 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 + + The zero element is never invertible:: + + sage: set_random_seed() + sage: J = random_eja().zero().inverse() + Traceback (most recent call last): + ... + ValueError: element is not invertible + + """ + 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:: + + sage: set_random_seed() + sage: J = random_eja() + sage: J.zero().is_invertible() + False + + """ + zero = self.parent().zero() + p = self.minimal_polynomial() + return not (p(zero) == zero) + + + 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:: + + sage: set_random_seed() + sage: random_eja().one().is_nilpotent() + 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:: + + sage: set_random_seed() + sage: J = random_eja() + sage: J.zero().is_regular() + False + + 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: n = ZZ.random_element(1,10) + sage: J = JordanSpinEJA(n) + 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:: + + sage: set_random_seed() + sage: J = random_eja() + sage: J.zero().degree() + 1 + sage: J.one().degree() + 1 + + Our implementation agrees with the definition:: + + sage: set_random_seed() + sage: x = random_eja().random_element() + sage: x.degree() == x.minimal_polynomial().degree() + True + + """ + return self.span_of_powers().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, + ....: random_eja) + + TESTS: + + The minimal polynomial of the identity and zero elements are + always the same:: + + sage: set_random_seed() + sage: J = random_eja() + 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:: + + sage: set_random_seed() + sage: n = ZZ.random_element(2,10) + sage: J = JordanSpinEJA(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: 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 + + """ + 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())) + return elt.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 + 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] + + :: + + sage: J = QuaternionHermitianEJA(3) + sage: J.one() + e0 + e9 + 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 = B[0].matrix_space() + return W.linear_combination(zip(self.vector(), B)) + + + 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 = J.random_element() + sage: y = J.random_element() + sage: x.operator()(y) == x*y + True + sage: y.operator()(x) == x*y + True + + """ + P = self.parent() + fda_elt = FiniteDimensionalAlgebraElement(P, self) + return FiniteDimensionalEuclideanJordanAlgebraOperator( + P, + P, + fda_elt.matrix().transpose() ) + + + 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: n = ZZ.random_element(1,10) + sage: J = JordanSpinEJA(n) + sage: x = J.random_element() + sage: x_vec = x.vector() + 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 = J.random_element() + sage: y = J.random_element() + 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 = QQ.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*x.operator()*Qex - Qx ) + True + + sage: 2*x.operator()*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 span_of_powers(self): + """ + Return the vector space spanned by successive powers of + this element. + """ + # The dimension of the subalgebra can't be greater than + # the big algebra, so just put everything into a list + # and let span() get rid of the excess. + # + # We do the extra ambient_vector_space() in case we're messing + # with polynomials and the direct parent is a module. + V = self.parent().vector_space() + return V.span( (self**d).vector() for d in xrange(V.dimension()) ) + + + def subalgebra_generated_by(self): + """ + Return the associative subalgebra of the parent EJA generated + by this element. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS:: + + sage: set_random_seed() + sage: x = random_eja().random_element() + sage: x.subalgebra_generated_by().is_associative() + True + + Squaring in the subalgebra should work the same as in + the superalgebra:: + + sage: set_random_seed() + sage: x = random_eja().random_element() + sage: u = x.subalgebra_generated_by().random_element() + sage: u.operator()(u) == u^2 + True + + """ + # First get the subspace spanned by the powers of myself... + V = self.span_of_powers() + F = self.base_ring() + + # Now figure out the entries of the right-multiplication + # matrix for the successive basis elements b0, b1,... of + # that subspace. + mats = [] + for b_right in V.basis(): + eja_b_right = self.parent()(b_right) + b_right_rows = [] + # The first row of the right-multiplication matrix by + # b1 is what we get if we apply that matrix to b1. The + # second row of the right multiplication matrix by b1 + # is what we get when we apply that matrix to b2... + # + # IMPORTANT: this assumes that all vectors are COLUMN + # vectors, unlike our superclass (which uses row vectors). + for b_left in V.basis(): + eja_b_left = self.parent()(b_left) + # Multiply in the original EJA, but then get the + # coordinates from the subalgebra in terms of its + # basis. + this_row = V.coordinates((eja_b_left*eja_b_right).vector()) + b_right_rows.append(this_row) + b_right_matrix = matrix(F, b_right_rows) + mats.append(b_right_matrix) + + # It's an algebra of polynomials in one element, and EJAs + # are power-associative. + # + # TODO: choose generator names intelligently. + # + # The rank is the highest possible degree of a minimal polynomial, + # and is bounded above by the dimension. We know in this case that + # there's an element whose minimal polynomial has the same degree + # as the space's dimension, so that must be its rank too. + return FiniteDimensionalEuclideanJordanAlgebra( + F, + mats, + V.dimension(), + assume_associative=True, + names='f') + + + 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:: + + sage: set_random_seed() + sage: J = random_eja() + 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.is_nilpotent(): + raise ValueError("this only works with non-nilpotent elements!") + + V = self.span_of_powers() + J = 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! + u = J(V.coordinates(self.vector())) + + # The image of the matrix of left-u^m-multiplication + # will be minimal for some natural number s... + s = 0 + minimal_dim = V.dimension() + for i in xrange(1, V.dimension()): + 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_coordinates = A.solve_right(u_next.vector()) + + # Now c_coordinates is the idempotent we want, but it's in + # the coordinate system of the subalgebra. + # + # We need the basis for J, but as elements of the parent algebra. + # + basis = [self.parent(v) for v in V.basis()] + return self.parent().linear_combination(zip(c_coordinates, basis)) + + + def trace(self): + """ + Return my trace, the sum of my eigenvalues. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + ....: RealCartesianProductEJA, + ....: random_eja) + + EXAMPLES:: + + 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 J.base_ring() + True + + """ + P = self.parent() + r = P.rank() + 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.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:: + + sage: set_random_seed() + sage: J = random_eja() + sage: x = J.random_element(); y = J.random_element() + sage: x.trace_inner_product(y) == y.trace_inner_product(x) + True + + The trace inner product is bilinear:: + + sage: set_random_seed() + sage: J = random_eja() + sage: x = J.random_element() + sage: y = J.random_element() + sage: z = J.random_element() + sage: a = QQ.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 + + The trace inner product satisfies the compatibility + condition in the definition of a Euclidean Jordan algebra:: + + sage: set_random_seed() + sage: J = random_eja() + sage: x = J.random_element() + sage: y = J.random_element() + sage: z = J.random_element() + 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() diff --git a/mjo/eja/eja_utils.py b/mjo/eja/eja_utils.py new file mode 100644 index 0000000..e8b7dc7 --- /dev/null +++ b/mjo/eja/eja_utils.py @@ -0,0 +1,9 @@ +from sage.matrix.constructor import matrix +from sage.modules.free_module_element import vector +from sage.functions.other import sqrt + +def _mat2vec(m): + return vector(m.base_ring(), m.list()) + +def _vec2mat(v): + return matrix(v.base_ring(), sqrt(v.degree()), v.list()) -- 2.44.2