X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_element.py;h=c98d9a2b64651562928a6489ef88e03fdd2b0dc9;hb=bc02bf48592e22d034310cfffef8fb2a062c0a43;hp=97c048dceb3e299e7a36ac1a15767ebb33af8fad;hpb=16dfa403c6eb709d3a5188a0f19919652b6a225d;p=sage.d.git diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py index 97c048d..c98d9a2 100644 --- a/mjo/eja/eja_element.py +++ b/mjo/eja/eja_element.py @@ -1,16 +1,12 @@ from sage.matrix.constructor import matrix +from sage.misc.cachefunc import cached_method 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_subalgebra', - 'FiniteDimensionalEuclideanJordanElementSubalgebra') -from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator -from mjo.eja.eja_utils import _mat2vec +from mjo.eja.eja_operator import FiniteDimensionalEJAOperator +from mjo.eja.eja_utils import _mat2vec, _scale -class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): +class FiniteDimensionalEJAElement(IndexedFreeModuleElement): """ An element of a Euclidean Jordan algebra. """ @@ -32,7 +28,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): Return ``self`` raised to the power ``n``. Jordan algebras are always power-associative; see for - example Faraut and Koranyi, Proposition II.1.2 (ii). + 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, @@ -78,7 +74,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): elif n == 1: return self else: - return (self.operator()**(n-1))(self) + return (self**(n-1))*self def apply_univariate_polynomial(self, p): @@ -94,7 +90,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): SETUP:: - sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA, + sage: from mjo.eja.eja_algebra import (HadamardEJA, ....: random_eja) EXAMPLES:: @@ -102,7 +98,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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 = HadamardEJA(5) sage: J.one().apply_univariate_polynomial(p) == 3*J.one() True @@ -111,7 +107,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): We should always get back an element of the algebra:: sage: set_random_seed() - sage: p = PolynomialRing(QQ, 't').random_element() + sage: p = PolynomialRing(AA, 't').random_element() sage: J = random_eja() sage: x = J.random_element() sage: x.apply_univariate_polynomial(p) in J @@ -135,7 +131,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): SETUP:: - sage: from mjo.eja.eja_algebra import RealCartesianProductEJA + sage: from mjo.eja.eja_algebra import HadamardEJA EXAMPLES: @@ -143,14 +139,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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 = HadamardEJA(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 = HadamardEJA(3) sage: J.zero().characteristic_polynomial() t^3 @@ -160,13 +156,28 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): to zero on that element:: sage: set_random_seed() - sage: x = RealCartesianProductEJA(3).random_element() + sage: x = HadamardEJA(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 = HadamardEJA(3) + sage: p1 = J.one().characteristic_polynomial() + sage: q1 = J.zero().characteristic_polynomial() + sage: b0,b1,b2 = J.gens() + sage: A = (b0 + 2*b1 + 3*b2).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() + p = self.parent().characteristic_polynomial_of() return p(*self.to_vector()) @@ -217,9 +228,9 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): Ditto for the quaternions:: - sage: J = QuaternionHermitianEJA(3) + sage: J = QuaternionHermitianEJA(2) sage: J.one().inner_product(J.one()) - 3 + 2 TESTS: @@ -228,9 +239,8 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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 + sage: x,y = J.random_elements(2) + sage: x.inner_product(y) in RLF True """ @@ -265,9 +275,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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: 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 @@ -277,9 +285,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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: x,y = random_eja().random_elements(2) sage: Lx = x.operator() sage: Ly = y.operator() sage: Lxx = (x*x).operator() @@ -291,10 +297,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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: x,y,z = random_eja().random_elements(3) sage: Lx = x.operator() sage: Ly = y.operator() sage: Lz = z.operator() @@ -308,10 +311,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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: u,y,z = random_eja().random_elements(3) sage: Lu = u.operator() sage: Ly = y.operator() sage: Lz = z.operator() @@ -340,12 +340,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + ....: TrivialEJA, + ....: RealSymmetricEJA, + ....: ComplexHermitianEJA, ....: random_eja) EXAMPLES:: sage: J = JordanSpinEJA(2) - sage: e0,e1 = J.gens() sage: x = sum( J.gens() ) sage: x.det() 0 @@ -353,11 +355,21 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): :: sage: J = JordanSpinEJA(3) - sage: e0,e1,e2 = J.gens() sage: x = sum( J.gens() ) sage: x.det() -1 + The determinant of the sole element in the rank-zero trivial + algebra is ``1``, by three paths of reasoning. First, its + characteristic polynomial is a constant ``1``, so the constant + term in that polynomial is ``1``. Second, the characteristic + polynomial evaluated at zero is again ``1``. And finally, the + (empty) product of its eigenvalues is likewise just unity:: + + sage: J = TrivialEJA() + sage: J.zero().det() + 1 + TESTS: An element is invertible if and only if its determinant is @@ -368,29 +380,67 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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 + + The determinant in real matrix algebras is the usual determinant:: + + sage: set_random_seed() + sage: X = matrix.random(QQ,3) + sage: X = X + X.T + sage: J1 = RealSymmetricEJA(3) + sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False) + sage: expected = X.det() + sage: actual1 = J1(X).det() + sage: actual2 = J2(X).det() + sage: actual1 == expected + True + sage: actual2 == expected + 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. + + if r == 0: + # Special case, since we don't get the a0=1 + # coefficient when the rank of the algebra + # is zero. + return P.base_ring().one() + + p = P._charpoly_coefficients()[0] + # The _charpoly_coeff function already adds the factor of -1 + # to ensure that _charpoly_coefficients()[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()) + @cached_method 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. + In general we appeal to the quadratic representation as in + Koecher's Theorem 12 in Chapter III, Section 5. But if the + parent algebra's "characteristic polynomial of" coefficients + happen to be cached, then we use Proposition II.2.4 in Faraut + and Korányi which gives a formula for the inverse based on the + characteristic polynomial and the Cayley-Hamilton theorem for + Euclidean Jordan algebras:: SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, + ....: JordanSpinEJA, ....: random_eja) EXAMPLES: @@ -399,20 +449,26 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): Example 11.11:: sage: set_random_seed() - sage: n = ZZ.random_element(1,10) - sage: J = JordanSpinEJA(n) + 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: x0 = x_vec[:1] 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: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar) + sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list()) + sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff 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): + ... + ZeroDivisionError: element is not invertible + TESTS: The identity element is its own inverse:: @@ -438,36 +494,63 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: (not x.is_invertible()) or (x.inverse().inverse() == x) True - The zero element is never invertible:: + 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().zero().inverse() - Traceback (most recent call last): - ... - ValueError: element is not invertible + sage: J = random_eja() + sage: x = J.random_element() + sage: (not x.operator().is_invertible()) or ( + ....: x.operator().inverse()(J.one()) == x.inverse() ) + True + Check that the fast (cached) and slow algorithms give the same + answer:: + + sage: set_random_seed() # long time + sage: J = random_eja(field=QQ, orthonormalize=False) # long time + sage: x = J.random_element() # long time + sage: while not x.is_invertible(): # long time + ....: x = J.random_element() # long time + sage: slow = x.inverse() # long time + sage: _ = J._charpoly_coefficients() # long time + sage: fast = x.inverse() # long time + sage: slow == fast # long time + True """ - if not self.is_invertible(): - raise ValueError("element is not invertible") - - return (~self.quadratic_representation())(self) - - + not_invertible_msg = "element is not invertible" + if self.parent()._charpoly_coefficients.is_in_cache(): + # We can invert using our charpoly if it will be fast to + # compute. If the coefficients are cached, our rank had + # better be too! + if self.det().is_zero(): + raise ZeroDivisionError(not_invertible_msg) + r = self.parent().rank() + a = self.characteristic_polynomial().coefficients(sparse=False) + return (-1)**(r+1)*sum(a[i+1]*self**i for i in range(r))/self.det() + + try: + inv = (~self.quadratic_representation())(self) + self.is_invertible.set_cache(True) + return inv + except ZeroDivisionError: + self.is_invertible.set_cache(False) + raise ZeroDivisionError(not_invertible_msg) + + + @cached_method 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. + If computing my determinant will be fast, we do so and compare + with zero (Proposition II.2.4 in Faraut and + Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi + reduces the problem to the invertibility of my quadratic + representation. SETUP:: @@ -482,19 +565,152 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: J.one().is_invertible() True - The zero element is never invertible:: + 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 + + Test that the fast (cached) and slow algorithms give the same + answer:: + + sage: set_random_seed() # long time + sage: J = random_eja(field=QQ, orthonormalize=False) # long time + sage: x = J.random_element() # long time + sage: slow = x.is_invertible() # long time + sage: _ = J._charpoly_coefficients() # long time + sage: fast = x.is_invertible() # long time + sage: slow == fast # long time + True + """ + if self.is_zero(): + if self.parent().is_trivial(): + return True + else: + return False + + if self.parent()._charpoly_coefficients.is_in_cache(): + # The determinant will be quicker than inverting the + # quadratic representation, most likely. + return (not self.det().is_zero()) + + # The easiest way to determine if I'm invertible is to try. + try: + inv = (~self.quadratic_representation())(self) + self.inverse.set_cache(inv) + return True + except ZeroDivisionError: + return False + + + 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) + 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) + sage: x = sum( i*J.monomial(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.zero().is_invertible() + 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 """ - # 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) + 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): @@ -524,10 +740,11 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): TESTS: - The identity element is never nilpotent:: + The identity element is never nilpotent, except in a trivial EJA:: sage: set_random_seed() - sage: random_eja().one().is_nilpotent() + sage: J = random_eja() + sage: J.one().is_nilpotent() and not J.is_trivial() False The additive identity is always nilpotent:: @@ -559,7 +776,9 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: J = JordanSpinEJA(5) sage: J.one().is_regular() False - sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity + sage: b0, b1, b2, b3, b4 = J.gens() + sage: b0 == J.one() + True sage: for x in J.gens(): ....: (J.one() + x).is_regular() False @@ -571,11 +790,11 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): TESTS: The zero element should never be regular, unless the parent - algebra has dimension one:: + 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() + sage: J.dimension() <= 1 or not J.zero().is_regular() True The unit element isn't regular unless the algebra happens to @@ -583,7 +802,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: set_random_seed() sage: J = random_eja() - sage: J.dimension() == 1 or not J.one().is_regular() + sage: J.dimension() <= 1 or not J.one().is_regular() True """ @@ -597,10 +816,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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:: @@ -612,30 +828,33 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: J = JordanSpinEJA(4) sage: J.one().degree() 1 - sage: e0,e1,e2,e3 = J.gens() - sage: (e0 - e1).degree() + sage: b0,b1,b2,b3 = J.gens() + sage: (b0 - b1).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: J = JordanSpinEJA.random_instance() + sage: n = J.dimension() sage: x = J.random_element() - sage: x == x.coefficient(0)*J.one() or x.degree() == 2 + sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one()) True TESTS: - The zero and unit elements are both of degree one:: + The zero and unit elements are both of degree one in nontrivial + algebras:: sage: set_random_seed() sage: J = random_eja() - sage: J.zero().degree() - 1 - sage: J.one().degree() - 1 + 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:: @@ -645,7 +864,59 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): True """ - return self.subalgebra_generated_by().dimension() + n = self.parent().dimension() + + if n == 0: + # The minimal polynomial is an empty product, i.e. the + # constant polynomial "1" having degree zero. + return 0 + elif self.is_zero(): + # The minimal polynomial of zero in a nontrivial algebra + # is "t", and is of degree one. + return 1 + elif n == 1: + # If this is a nonzero element of a nontrivial algebra, it + # has degree at least one. It follows that, in an algebra + # of dimension one, the degree must be actually one. + return 1 + + # BEWARE: The subalgebra_generated_by() method uses the result + # of this method to construct a basis for the subalgebra. That + # means, in particular, that we cannot implement this method + # as ``self.subalgebra_generated_by().dimension()``. + + # Algorithm: keep appending (vector representations of) powers + # self as rows to a matrix and echelonizing it. When its rank + # stops increasing, we've reached a redundancy. + + # Given the special cases above, we can assume that "self" is + # nonzero, the algebra is nontrivial, and that its dimension + # is at least two. + M = matrix([(self.parent().one()).to_vector()]) + old_rank = 1 + + # Specifying the row-reduction algorithm can e.g. help over + # AA because it avoids the RecursionError that gets thrown + # when we have to look too hard for a root. + # + # Beware: QQ supports an entirely different set of "algorithm" + # keywords than do AA and RR. + algo = None + from sage.rings.all import QQ + if self.parent().base_ring() is not QQ: + algo = "scaled_partial_pivoting" + + for d in range(1,n): + M = matrix(M.rows() + [(self**d).to_vector()]) + M.echelonize(algo) + new_rank = M.rank() + if new_rank == old_rank: + return new_rank + else: + old_rank = new_rank + + return n + def left_matrix(self): @@ -673,18 +944,38 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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:: + always the same, except in trivial algebras where the minimal + polynomial of the unit/zero element is ``1``:: sage: set_random_seed() sage: J = random_eja() - sage: J.one().minimal_polynomial() + sage: mu = J.one().minimal_polynomial() + sage: t = mu.parent().gen() + sage: mu + int(J.is_trivial())*(t-2) t - 1 - sage: J.zero().minimal_polynomial() + sage: mu = J.zero().minimal_polynomial() + sage: t = mu.parent().gen() + sage: mu + int(J.is_trivial())*(t-1) t The degree of an element is (by one definition) the degree @@ -698,10 +989,12 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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:: + 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 = ZZ.random_element(2,10) + sage: n_max = max(2, JordanSpinEJA._max_random_instance_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(): @@ -722,63 +1015,150 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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_random_instance_size() + sage: n = ZZ.random_element(1, n_max) + sage: J1 = RealSymmetricEJA(n) + sage: J2 = RealSymmetricEJA(n,orthonormalize=False) + sage: X = random_matrix(AA,n) + sage: X = X*X.transpose() + sage: x1 = J1(X) + sage: x2 = J2(X) + sage: x1.minimal_polynomial() == x2.minimal_polynomial() + True + """ - A = self.subalgebra_generated_by() - return A(self).operator().minimal_polynomial() + if self.is_zero(): + # 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() + + # If we don't orthonormalize the subalgebra's basis, then the + # first two monomials in the subalgebra will be self^0 and + # self^1... assuming that self^1 is not a scalar multiple of + # self^0 (the unit element). We special case these to avoid + # having to solve a system to coerce self into the subalgebra. + A = self.subalgebra_generated_by(orthonormalize=False) + + if A.dimension() == 1: + # Does a solve to find the scalar multiple alpha such that + # alpha*unit = self. We have to do this because the basis + # for the subalgebra will be [ self^0 ], and not [ self^1 ]! + unit = self.parent().one() + alpha = self.to_vector() / unit.to_vector() + return (unit.operator()*alpha).minimal_polynomial() + else: + # If the dimension of the subalgebra is >= 2, then we just + # use the second basis element. + return A.monomial(1).operator().minimal_polynomial() - def natural_representation(self): + def to_matrix(self): """ - Return a more-natural representation of this element. + Return an (often more natural) representation of this element as a + matrix. - 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. + Every finite-dimensional Euclidean Jordan Algebra is a direct + sum of five simple algebras, four of which comprise Hermitian + matrices. This method returns a "natural" matrix + representation of this element as either a Hermitian matrix or + column vector. SETUP:: sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, - ....: QuaternionHermitianEJA) + ....: HadamardEJA, + ....: QuaternionHermitianEJA, + ....: RealSymmetricEJA) 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] + b0 + b3 + b8 + sage: J.one().to_matrix() + +---+---+---+ + | 1 | 0 | 0 | + +---+---+---+ + | 0 | 1 | 0 | + +---+---+---+ + | 0 | 0 | 1 | + +---+---+---+ :: - sage: J = QuaternionHermitianEJA(3) + sage: J = QuaternionHermitianEJA(2) 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] + b0 + b5 + sage: J.one().to_matrix() + +---+---+ + | 1 | 0 | + +---+---+ + | 0 | 1 | + +---+---+ + + This also works in Cartesian product algebras:: + + sage: J1 = HadamardEJA(1) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: x = sum(J.gens()) + sage: x.to_matrix()[0] + [1] + sage: x.to_matrix()[1] + [ 1 0.7071067811865475?] + [0.7071067811865475? 1] """ - B = self.parent().natural_basis() - W = B[0].matrix_space() - return W.linear_combination(zip(B,self.to_vector())) + B = self.parent().matrix_basis() + W = self.parent().matrix_space() + + if hasattr(W, 'cartesian_factors'): + # Aaaaand linear combinations don't work in Cartesian + # product spaces, even though they provide a method with + # that name. This is hidden behind an "if" because the + # _scale() function is slow. + pairs = zip(B, self.to_vector()) + return W.sum( _scale(b, alpha) for (b,alpha) in pairs ) + else: + # This is just a manual "from_vector()", but of course + # matrix spaces aren't vector spaces in sage, so they + # don't have a from_vector() method. + return W.linear_combination( zip(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, + ....: HadamardEJA) + + EXAMPLES:: + + sage: J = HadamardEJA(2) + sage: x = sum(J.gens()) + sage: x.norm() + 1.414213562373095? + + :: + + sage: J = JordanSpinEJA(4) + sage: x = sum(J.gens()) + sage: x.norm() + 2 + + """ + return self.inner_product(self).sqrt() def operator(self): @@ -794,8 +1174,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: set_random_seed() sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() + sage: x,y = J.random_elements(2) sage: x.operator()(y) == x*y True sage: y.operator()(x) == x*y @@ -805,10 +1184,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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() ) + return FiniteDimensionalEJAOperator(P, P, L.matrix() ) def quadratic_representation(self, other=None): @@ -826,19 +1202,20 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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 = JordanSpinEJA.random_instance().random_element() sage: x_vec = x.to_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 = matrix.identity(x.base_ring(), 0) + sage: n = x_vec.degree() + sage: if n > 0: + ....: x0 = x_vec[0] + ....: x_bar = x_vec[1:] + ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)]) + ....: B = 2*x0*x_bar.row() + ....: C = 2*x0*x_bar.column() + ....: D = matrix.identity(x.base_ring(), n-1) + ....: D = (x0^2 - x_bar.inner_product(x_bar))*D + ....: D = D + 2*x_bar.tensor_product(x_bar) + ....: Q = matrix.block(2,2,[A,B,C,D]) sage: Q == x.quadratic_representation().matrix() True @@ -846,8 +1223,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: set_random_seed() sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() + sage: x,y = J.random_elements(2) sage: Lx = x.operator() sage: Lxx = (x*x).operator() sage: Qx = x.quadratic_representation() @@ -864,7 +1240,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): Property 2 (multiply on the right for :trac:`28272`): - sage: alpha = QQ.random_element() + sage: alpha = J.base_ring().random_element() sage: (alpha*x).quadratic_representation() == Qx*(alpha^2) True @@ -892,10 +1268,10 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: not x.is_invertible() or ( ....: x.quadratic_representation(x.inverse())*Qx ....: == - ....: 2*x.operator()*Qex - Qx ) + ....: 2*Lx*Qex - Qx ) True - sage: 2*x.operator()*Qex - Qx == Lxx + sage: 2*Lx*Qex - Qx == Lxx True Property 5: @@ -932,15 +1308,108 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): + 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) + sage: J.one() + b0 + b2 + b5 + sage: J.one().spectral_decomposition() + [(1, b0 + b2 + b5)] + sage: J.zero().spectral_decomposition() + [(0, b0 + b2 + b5)] + + TESTS:: + + sage: J = RealSymmetricEJA(4) + 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 + + The spectral decomposition should work in subalgebras, too:: + + sage: J = RealSymmetricEJA(4) + sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens() + sage: A = 2*b5 - 2*b8 + sage: (lambda1, c1) = A.spectral_decomposition()[1] + sage: (J0, J5, J1) = J.peirce_decomposition(c1) + sage: (f0, f1, f2) = J1.gens() + sage: f0.spectral_decomposition() + [(0, c2), (1, c0)] + + """ + A = self.subalgebra_generated_by(orthonormalize=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): + def subalgebra_generated_by(self, **kwargs): """ 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 + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, + ....: RealSymmetricEJA) + + EXAMPLES: + + We can create subalgebras of Cartesian product EJAs that are not + themselves Cartesian product EJAs (they're just "regular" EJAs):: + + sage: J1 = HadamardEJA(3) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J.one().subalgebra_generated_by() + Euclidean Jordan algebra of dimension 1 over Algebraic Real Field TESTS: @@ -949,9 +1418,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: set_random_seed() sage: x0 = random_eja().random_element() sage: A = x0.subalgebra_generated_by() - sage: x = A.random_element() - sage: y = A.random_element() - sage: z = A.random_element() + sage: x,y,z = A.random_elements(3) sage: (x*y)*z == x*(y*z) True @@ -964,17 +1431,25 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: A(x^2) == A(x)*A(x) True - The subalgebra generated by the zero element is trivial:: + 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 - Euclidean Jordan algebra of dimension 0 over Rational Field - sage: A.one() - 0 + sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1 + True """ - return FiniteDimensionalEuclideanJordanElementSubalgebra(self) + powers = tuple( self**k for k in range(self.degree()) ) + A = self.parent().subalgebra(powers, + associative=True, + check_field=False, + check_axioms=False, + **kwargs) + A.one.set_cache(A(self.parent().one())) + return A def subalgebra_idempotent(self): @@ -986,18 +1461,25 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: from mjo.eja.eja_algebra import random_eja - TESTS:: + TESTS: + + Ensure that we can find an idempotent in a non-trivial algebra + where there are non-nilpotent elements, or that we get the dumb + solution in the trivial algebra:: sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() - sage: while x.is_nilpotent(): + sage: while x.is_nilpotent() and not J.is_trivial(): ....: 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!") @@ -1008,7 +1490,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): # will be minimal for some natural number s... s = 0 minimal_dim = J.dimension() - for i in xrange(1, minimal_dim): + for i in range(1, minimal_dim): this_dim = (u**i).operator().matrix().image().dimension() if this_dim < minimal_dim: minimal_dim = this_dim @@ -1037,14 +1519,23 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): """ 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, + ....: HadamardEJA, + ....: TrivialEJA, ....: random_eja) EXAMPLES:: + sage: J = TrivialEJA() + sage: J.zero().trace() + 0 + + :: sage: J = JordanSpinEJA(3) sage: x = sum(J.gens()) sage: x.trace() @@ -1052,7 +1543,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): :: - sage: J = RealCartesianProductEJA(5) + sage: J = HadamardEJA(5) sage: J.one().trace() 5 @@ -1062,13 +1553,28 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: set_random_seed() sage: J = random_eja() - sage: J.random_element().trace() in J.base_ring() + sage: J.random_element().trace() in RLF + True + + The trace is linear:: + + sage: set_random_seed() + sage: J = random_eja() + sage: x,y = J.random_elements(2) + sage: alpha = J.base_ring().random_element() + sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace() True """ P = self.parent() r = P.rank() - p = P._charpoly_coeff(r-1) + + if r == 0: + # Special case for the trivial algebra where + # the trace is an empty sum. + return P.base_ring().zero() + + p = P._charpoly_coefficients()[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, @@ -1086,22 +1592,16 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): TESTS: - The trace inner product is commutative:: + The trace inner product is commutative, bilinear, and associative:: sage: set_random_seed() sage: J = random_eja() - sage: x = J.random_element(); y = J.random_element() + sage: x,y,z = J.random_elements(3) + sage: # commutative 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: # 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) ) @@ -1112,15 +1612,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): ....: 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: # associative sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z) True @@ -1129,3 +1621,30 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): 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, + ....: HadamardEJA) + + EXAMPLES:: + + sage: J = HadamardEJA(2) + sage: x = sum(J.gens()) + sage: x.trace_norm() + 1.414213562373095? + + :: + + sage: J = JordanSpinEJA(4) + sage: x = sum(J.gens()) + sage: x.trace_norm() + 2.828427124746190? + + """ + return self.trace_inner_product(self).sqrt()