1 from sage
.matrix
.constructor
import matrix
2 from sage
.modules
.free_module
import VectorSpace
3 from sage
.modules
.with_basis
.indexed_element
import IndexedFreeModuleElement
5 from mjo
.eja
.eja_operator
import FiniteDimensionalEuclideanJordanAlgebraOperator
6 from mjo
.eja
.eja_utils
import _mat2vec
8 class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement
):
10 An element of a Euclidean Jordan algebra.
15 Oh man, I should not be doing this. This hides the "disabled"
16 methods ``left_matrix`` and ``matrix`` from introspection;
17 in particular it removes them from tab-completion.
19 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
27 Return ``self`` raised to the power ``n``.
29 Jordan algebras are always power-associative; see for
30 example Faraut and Korányi, Proposition II.1.2 (ii).
32 We have to override this because our superclass uses row
33 vectors instead of column vectors! We, on the other hand,
34 assume column vectors everywhere.
38 sage: from mjo.eja.eja_algebra import random_eja
42 The definition of `x^2` is the unambiguous `x*x`::
44 sage: set_random_seed()
45 sage: x = random_eja().random_element()
49 A few examples of power-associativity::
51 sage: set_random_seed()
52 sage: x = random_eja().random_element()
53 sage: x*(x*x)*(x*x) == x^5
55 sage: (x*x)*(x*x*x) == x^5
58 We also know that powers operator-commute (Koecher, Chapter
61 sage: set_random_seed()
62 sage: x = random_eja().random_element()
63 sage: m = ZZ.random_element(0,10)
64 sage: n = ZZ.random_element(0,10)
65 sage: Lxm = (x^m).operator()
66 sage: Lxn = (x^n).operator()
67 sage: Lxm*Lxn == Lxn*Lxm
72 return self
.parent().one()
76 return (self
**(n
-1))*self
79 def apply_univariate_polynomial(self
, p
):
81 Apply the univariate polynomial ``p`` to this element.
83 A priori, SageMath won't allow us to apply a univariate
84 polynomial to an element of an EJA, because we don't know
85 that EJAs are rings (they are usually not associative). Of
86 course, we know that EJAs are power-associative, so the
87 operation is ultimately kosher. This function sidesteps
88 the CAS to get the answer we want and expect.
92 sage: from mjo.eja.eja_algebra import (HadamardEJA,
97 sage: R = PolynomialRing(QQ, 't')
99 sage: p = t^4 - t^3 + 5*t - 2
100 sage: J = HadamardEJA(5)
101 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
106 We should always get back an element of the algebra::
108 sage: set_random_seed()
109 sage: p = PolynomialRing(AA, 't').random_element()
110 sage: J = random_eja()
111 sage: x = J.random_element()
112 sage: x.apply_univariate_polynomial(p) in J
116 if len(p
.variables()) > 1:
117 raise ValueError("not a univariate polynomial")
120 # Convert the coeficcients to the parent's base ring,
121 # because a priori they might live in an (unnecessarily)
122 # larger ring for which P.sum() would fail below.
123 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
124 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
127 def characteristic_polynomial(self
):
129 Return the characteristic polynomial of this element.
133 sage: from mjo.eja.eja_algebra import HadamardEJA
137 The rank of `R^3` is three, and the minimal polynomial of
138 the identity element is `(t-1)` from which it follows that
139 the characteristic polynomial should be `(t-1)^3`::
141 sage: J = HadamardEJA(3)
142 sage: J.one().characteristic_polynomial()
143 t^3 - 3*t^2 + 3*t - 1
145 Likewise, the characteristic of the zero element in the
146 rank-three algebra `R^{n}` should be `t^{3}`::
148 sage: J = HadamardEJA(3)
149 sage: J.zero().characteristic_polynomial()
154 The characteristic polynomial of an element should evaluate
155 to zero on that element::
157 sage: set_random_seed()
158 sage: x = HadamardEJA(3).random_element()
159 sage: p = x.characteristic_polynomial()
160 sage: x.apply_univariate_polynomial(p)
163 The characteristic polynomials of the zero and unit elements
164 should be what we think they are in a subalgebra, too::
166 sage: J = HadamardEJA(3)
167 sage: p1 = J.one().characteristic_polynomial()
168 sage: q1 = J.zero().characteristic_polynomial()
169 sage: e0,e1,e2 = J.gens()
170 sage: A = (e0 + 2*e1 + 3*e2).subalgebra_generated_by() # dim 3
171 sage: p2 = A.one().characteristic_polynomial()
172 sage: q2 = A.zero().characteristic_polynomial()
179 p
= self
.parent().characteristic_polynomial_of()
180 return p(*self
.to_vector())
183 def inner_product(self
, other
):
185 Return the parent algebra's inner product of myself and ``other``.
189 sage: from mjo.eja.eja_algebra import (
190 ....: ComplexHermitianEJA,
192 ....: QuaternionHermitianEJA,
193 ....: RealSymmetricEJA,
198 The inner product in the Jordan spin algebra is the usual
199 inner product on `R^n` (this example only works because the
200 basis for the Jordan algebra is the standard basis in `R^n`)::
202 sage: J = JordanSpinEJA(3)
203 sage: x = vector(QQ,[1,2,3])
204 sage: y = vector(QQ,[4,5,6])
205 sage: x.inner_product(y)
207 sage: J.from_vector(x).inner_product(J.from_vector(y))
210 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
211 multiplication is the usual matrix multiplication in `S^n`,
212 so the inner product of the identity matrix with itself
215 sage: J = RealSymmetricEJA(3)
216 sage: J.one().inner_product(J.one())
219 Likewise, the inner product on `C^n` is `<X,Y> =
220 Re(trace(X*Y))`, where we must necessarily take the real
221 part because the product of Hermitian matrices may not be
224 sage: J = ComplexHermitianEJA(3)
225 sage: J.one().inner_product(J.one())
228 Ditto for the quaternions::
230 sage: J = QuaternionHermitianEJA(3)
231 sage: J.one().inner_product(J.one())
236 Ensure that we can always compute an inner product, and that
237 it gives us back a real number::
239 sage: set_random_seed()
240 sage: J = random_eja()
241 sage: x,y = J.random_elements(2)
242 sage: x.inner_product(y) in RLF
248 raise TypeError("'other' must live in the same algebra")
250 return P
.inner_product(self
, other
)
253 def operator_commutes_with(self
, other
):
255 Return whether or not this element operator-commutes
260 sage: from mjo.eja.eja_algebra import random_eja
264 The definition of a Jordan algebra says that any element
265 operator-commutes with its square::
267 sage: set_random_seed()
268 sage: x = random_eja().random_element()
269 sage: x.operator_commutes_with(x^2)
274 Test Lemma 1 from Chapter III of Koecher::
276 sage: set_random_seed()
277 sage: u,v = random_eja().random_elements(2)
278 sage: lhs = u.operator_commutes_with(u*v)
279 sage: rhs = v.operator_commutes_with(u^2)
283 Test the first polarization identity from my notes, Koecher
284 Chapter III, or from Baes (2.3)::
286 sage: set_random_seed()
287 sage: x,y = random_eja().random_elements(2)
288 sage: Lx = x.operator()
289 sage: Ly = y.operator()
290 sage: Lxx = (x*x).operator()
291 sage: Lxy = (x*y).operator()
292 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
295 Test the second polarization identity from my notes or from
298 sage: set_random_seed()
299 sage: x,y,z = random_eja().random_elements(3)
300 sage: Lx = x.operator()
301 sage: Ly = y.operator()
302 sage: Lz = z.operator()
303 sage: Lzy = (z*y).operator()
304 sage: Lxy = (x*y).operator()
305 sage: Lxz = (x*z).operator()
306 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
309 Test the third polarization identity from my notes or from
312 sage: set_random_seed()
313 sage: u,y,z = random_eja().random_elements(3)
314 sage: Lu = u.operator()
315 sage: Ly = y.operator()
316 sage: Lz = z.operator()
317 sage: Lzy = (z*y).operator()
318 sage: Luy = (u*y).operator()
319 sage: Luz = (u*z).operator()
320 sage: Luyz = (u*(y*z)).operator()
321 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
322 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
323 sage: bool(lhs == rhs)
327 if not other
in self
.parent():
328 raise TypeError("'other' must live in the same algebra")
337 Return my determinant, the product of my eigenvalues.
341 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
343 ....: RealSymmetricEJA,
344 ....: ComplexHermitianEJA,
349 sage: J = JordanSpinEJA(2)
350 sage: e0,e1 = J.gens()
351 sage: x = sum( J.gens() )
357 sage: J = JordanSpinEJA(3)
358 sage: e0,e1,e2 = J.gens()
359 sage: x = sum( J.gens() )
363 The determinant of the sole element in the rank-zero trivial
364 algebra is ``1``, by three paths of reasoning. First, its
365 characteristic polynomial is a constant ``1``, so the constant
366 term in that polynomial is ``1``. Second, the characteristic
367 polynomial evaluated at zero is again ``1``. And finally, the
368 (empty) product of its eigenvalues is likewise just unity::
370 sage: J = TrivialEJA()
376 An element is invertible if and only if its determinant is
379 sage: set_random_seed()
380 sage: x = random_eja().random_element()
381 sage: x.is_invertible() == (x.det() != 0)
384 Ensure that the determinant is multiplicative on an associative
385 subalgebra as in Faraut and Korányi's Proposition II.2.2::
387 sage: set_random_seed()
388 sage: J = random_eja().random_element().subalgebra_generated_by()
389 sage: x,y = J.random_elements(2)
390 sage: (x*y).det() == x.det()*y.det()
393 The determinant in matrix algebras is just the usual determinant::
395 sage: set_random_seed()
396 sage: X = matrix.random(QQ,3)
398 sage: J1 = RealSymmetricEJA(3)
399 sage: J2 = RealSymmetricEJA(3,QQ,orthonormalize=False)
400 sage: expected = X.det()
401 sage: actual1 = J1(X).det()
402 sage: actual2 = J2(X).det()
403 sage: actual1 == expected
405 sage: actual2 == expected
410 sage: set_random_seed()
411 sage: J1 = ComplexHermitianEJA(3)
412 sage: J2 = ComplexHermitianEJA(3,field=QQ,orthonormalize=False)
413 sage: X = matrix.random(GaussianIntegers(),3)
415 sage: expected = AA(X.det())
416 sage: actual1 = J1(J1.real_embed(X)).det()
417 sage: actual2 = J2(J2.real_embed(X)).det()
418 sage: expected == actual1
420 sage: expected == actual2
428 # Special case, since we don't get the a0=1
429 # coefficient when the rank of the algebra
431 return P
.base_ring().one()
433 p
= P
._charpoly
_coefficients
()[0]
434 # The _charpoly_coeff function already adds the factor of -1
435 # to ensure that _charpoly_coefficients()[0] is really what
436 # appears in front of t^{0} in the charpoly. However, we want
437 # (-1)^r times THAT for the determinant.
438 return ((-1)**r
)*p(*self
.to_vector())
443 Return the Jordan-multiplicative inverse of this element.
447 We appeal to the quadratic representation as in Koecher's
448 Theorem 12 in Chapter III, Section 5.
452 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
458 The inverse in the spin factor algebra is given in Alizadeh's
461 sage: set_random_seed()
462 sage: J = JordanSpinEJA.random_instance()
463 sage: x = J.random_element()
464 sage: while not x.is_invertible():
465 ....: x = J.random_element()
466 sage: x_vec = x.to_vector()
468 sage: x_bar = x_vec[1:]
469 sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
470 sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
471 sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
472 sage: x.inverse() == J.from_vector(x_inverse)
475 Trying to invert a non-invertible element throws an error:
477 sage: JordanSpinEJA(3).zero().inverse()
478 Traceback (most recent call last):
480 ValueError: element is not invertible
484 The identity element is its own inverse::
486 sage: set_random_seed()
487 sage: J = random_eja()
488 sage: J.one().inverse() == J.one()
491 If an element has an inverse, it acts like one::
493 sage: set_random_seed()
494 sage: J = random_eja()
495 sage: x = J.random_element()
496 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
499 The inverse of the inverse is what we started with::
501 sage: set_random_seed()
502 sage: J = random_eja()
503 sage: x = J.random_element()
504 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
507 Proposition II.2.3 in Faraut and Korányi says that the inverse
508 of an element is the inverse of its left-multiplication operator
509 applied to the algebra's identity, when that inverse exists::
511 sage: set_random_seed()
512 sage: J = random_eja()
513 sage: x = J.random_element()
514 sage: (not x.operator().is_invertible()) or (
515 ....: x.operator().inverse()(J.one()) == x.inverse() )
518 Proposition II.2.4 in Faraut and Korányi gives a formula for
519 the inverse based on the characteristic polynomial and the
520 Cayley-Hamilton theorem for Euclidean Jordan algebras::
522 sage: set_random_seed()
523 sage: J = ComplexHermitianEJA(3)
524 sage: x = J.random_element()
525 sage: while not x.is_invertible():
526 ....: x = J.random_element()
528 sage: a = x.characteristic_polynomial().coefficients(sparse=False)
529 sage: expected = (-1)^(r+1)/x.det()
530 sage: expected *= sum( a[i+1]*x^i for i in range(r) )
531 sage: x.inverse() == expected
535 if not self
.is_invertible():
536 raise ValueError("element is not invertible")
538 if self
.parent()._charpoly
_coefficients
.is_in_cache():
539 # We can invert using our charpoly if it will be fast to
540 # compute. If the coefficients are cached, our rank had
542 r
= self
.parent().rank()
543 a
= self
.characteristic_polynomial().coefficients(sparse
=False)
544 return (-1)**(r
+1)*sum(a
[i
+1]*self
**i
for i
in range(r
))/self
.det()
546 return (~self
.quadratic_representation())(self
)
549 def is_invertible(self
):
551 Return whether or not this element is invertible.
555 The usual way to do this is to check if the determinant is
556 zero, but we need the characteristic polynomial for the
557 determinant. The minimal polynomial is a lot easier to get,
558 so we use Corollary 2 in Chapter V of Koecher to check
559 whether or not the parent algebra's zero element is a root
560 of this element's minimal polynomial.
562 That is... unless the coefficients of our algebra's
563 "characteristic polynomial of" function are already cached!
564 In that case, we just use the determinant (which will be fast
567 Beware that we can't use the superclass method, because it
568 relies on the algebra being associative.
572 sage: from mjo.eja.eja_algebra import random_eja
576 The identity element is always invertible::
578 sage: set_random_seed()
579 sage: J = random_eja()
580 sage: J.one().is_invertible()
583 The zero element is never invertible in a non-trivial algebra::
585 sage: set_random_seed()
586 sage: J = random_eja()
587 sage: (not J.is_trivial()) and J.zero().is_invertible()
592 if self
.parent().is_trivial():
597 if self
.parent()._charpoly
_coefficients
.is_in_cache():
598 # The determinant will be quicker than computing the minimal
599 # polynomial from scratch, most likely.
600 return (not self
.det().is_zero())
602 # In fact, we only need to know if the constant term is non-zero,
603 # so we can pass in the field's zero element instead.
604 zero
= self
.base_ring().zero()
605 p
= self
.minimal_polynomial()
606 return not (p(zero
) == zero
)
609 def is_primitive_idempotent(self
):
611 Return whether or not this element is a primitive (or minimal)
614 A primitive idempotent is a non-zero idempotent that is not
615 the sum of two other non-zero idempotents. Remark 2.7.15 in
616 Baes shows that this is what he refers to as a "minimal
619 An element of a Euclidean Jordan algebra is a minimal idempotent
620 if it :meth:`is_idempotent` and if its Peirce subalgebra
621 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
626 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
627 ....: RealSymmetricEJA,
633 This method is sloooooow.
637 The spectral decomposition of a non-regular element should always
638 contain at least one non-minimal idempotent::
640 sage: J = RealSymmetricEJA(3)
641 sage: x = sum(J.gens())
644 sage: [ c.is_primitive_idempotent()
645 ....: for (l,c) in x.spectral_decomposition() ]
648 On the other hand, the spectral decomposition of a regular
649 element should always be in terms of minimal idempotents::
651 sage: J = JordanSpinEJA(4)
652 sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) )
655 sage: [ c.is_primitive_idempotent()
656 ....: for (l,c) in x.spectral_decomposition() ]
661 The identity element is minimal only in an EJA of rank one::
663 sage: set_random_seed()
664 sage: J = random_eja()
665 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
668 A non-idempotent cannot be a minimal idempotent::
670 sage: set_random_seed()
671 sage: J = JordanSpinEJA(4)
672 sage: x = J.random_element()
673 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
676 Proposition 2.7.19 in Baes says that an element is a minimal
677 idempotent if and only if it's idempotent with trace equal to
680 sage: set_random_seed()
681 sage: J = JordanSpinEJA(4)
682 sage: x = J.random_element()
683 sage: expected = (x.is_idempotent() and x.trace() == 1)
684 sage: actual = x.is_primitive_idempotent()
685 sage: actual == expected
688 Primitive idempotents must be non-zero::
690 sage: set_random_seed()
691 sage: J = random_eja()
692 sage: J.zero().is_idempotent()
694 sage: J.zero().is_primitive_idempotent()
697 As a consequence of the fact that primitive idempotents must
698 be non-zero, there are no primitive idempotents in a trivial
699 Euclidean Jordan algebra::
701 sage: J = TrivialEJA()
702 sage: J.one().is_idempotent()
704 sage: J.one().is_primitive_idempotent()
708 if not self
.is_idempotent():
714 (_
,_
,J1
) = self
.parent().peirce_decomposition(self
)
715 return (J1
.dimension() == 1)
718 def is_nilpotent(self
):
720 Return whether or not some power of this element is zero.
724 We use Theorem 5 in Chapter III of Koecher, which says that
725 an element ``x`` is nilpotent if and only if ``x.operator()``
726 is nilpotent. And it is a basic fact of linear algebra that
727 an operator on an `n`-dimensional space is nilpotent if and
728 only if, when raised to the `n`th power, it equals the zero
729 operator (for example, see Axler Corollary 8.8).
733 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
738 sage: J = JordanSpinEJA(3)
739 sage: x = sum(J.gens())
740 sage: x.is_nilpotent()
745 The identity element is never nilpotent, except in a trivial EJA::
747 sage: set_random_seed()
748 sage: J = random_eja()
749 sage: J.one().is_nilpotent() and not J.is_trivial()
752 The additive identity is always nilpotent::
754 sage: set_random_seed()
755 sage: random_eja().zero().is_nilpotent()
760 zero_operator
= P
.zero().operator()
761 return self
.operator()**P
.dimension() == zero_operator
764 def is_regular(self
):
766 Return whether or not this is a regular element.
770 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
775 The identity element always has degree one, but any element
776 linearly-independent from it is regular::
778 sage: J = JordanSpinEJA(5)
779 sage: J.one().is_regular()
781 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
782 sage: for x in J.gens():
783 ....: (J.one() + x).is_regular()
792 The zero element should never be regular, unless the parent
793 algebra has dimension less than or equal to one::
795 sage: set_random_seed()
796 sage: J = random_eja()
797 sage: J.dimension() <= 1 or not J.zero().is_regular()
800 The unit element isn't regular unless the algebra happens to
801 consist of only its scalar multiples::
803 sage: set_random_seed()
804 sage: J = random_eja()
805 sage: J.dimension() <= 1 or not J.one().is_regular()
809 return self
.degree() == self
.parent().rank()
814 Return the degree of this element, which is defined to be
815 the degree of its minimal polynomial.
819 For now, we skip the messy minimal polynomial computation
820 and instead return the dimension of the vector space spanned
821 by the powers of this element. The latter is a bit more
822 straightforward to compute.
826 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
831 sage: J = JordanSpinEJA(4)
832 sage: J.one().degree()
834 sage: e0,e1,e2,e3 = J.gens()
835 sage: (e0 - e1).degree()
838 In the spin factor algebra (of rank two), all elements that
839 aren't multiples of the identity are regular::
841 sage: set_random_seed()
842 sage: J = JordanSpinEJA.random_instance()
843 sage: n = J.dimension()
844 sage: x = J.random_element()
845 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
850 The zero and unit elements are both of degree one in nontrivial
853 sage: set_random_seed()
854 sage: J = random_eja()
855 sage: d = J.zero().degree()
856 sage: (J.is_trivial() and d == 0) or d == 1
858 sage: d = J.one().degree()
859 sage: (J.is_trivial() and d == 0) or d == 1
862 Our implementation agrees with the definition::
864 sage: set_random_seed()
865 sage: x = random_eja().random_element()
866 sage: x.degree() == x.minimal_polynomial().degree()
870 if self
.is_zero() and not self
.parent().is_trivial():
871 # The minimal polynomial of zero in a nontrivial algebra
872 # is "t"; in a trivial algebra it's "1" by convention
873 # (it's an empty product).
875 return self
.subalgebra_generated_by().dimension()
878 def left_matrix(self
):
880 Our parent class defines ``left_matrix`` and ``matrix``
881 methods whose names are misleading. We don't want them.
883 raise NotImplementedError("use operator().matrix() instead")
888 def minimal_polynomial(self
):
890 Return the minimal polynomial of this element,
891 as a function of the variable `t`.
895 We restrict ourselves to the associative subalgebra
896 generated by this element, and then return the minimal
897 polynomial of this element's operator matrix (in that
898 subalgebra). This works by Baes Proposition 2.3.16.
902 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
903 ....: RealSymmetricEJA,
909 Keeping in mind that the polynomial ``1`` evaluates the identity
910 element (also the zero element) of the trivial algebra, it is clear
911 that the polynomial ``1`` is the minimal polynomial of the only
912 element in a trivial algebra::
914 sage: J = TrivialEJA()
915 sage: J.one().minimal_polynomial()
917 sage: J.zero().minimal_polynomial()
922 The minimal polynomial of the identity and zero elements are
923 always the same, except in trivial algebras where the minimal
924 polynomial of the unit/zero element is ``1``::
926 sage: set_random_seed()
927 sage: J = random_eja()
928 sage: mu = J.one().minimal_polynomial()
929 sage: t = mu.parent().gen()
930 sage: mu + int(J.is_trivial())*(t-2)
932 sage: mu = J.zero().minimal_polynomial()
933 sage: t = mu.parent().gen()
934 sage: mu + int(J.is_trivial())*(t-1)
937 The degree of an element is (by one definition) the degree
938 of its minimal polynomial::
940 sage: set_random_seed()
941 sage: x = random_eja().random_element()
942 sage: x.degree() == x.minimal_polynomial().degree()
945 The minimal polynomial and the characteristic polynomial coincide
946 and are known (see Alizadeh, Example 11.11) for all elements of
947 the spin factor algebra that aren't scalar multiples of the
948 identity. We require the dimension of the algebra to be at least
949 two here so that said elements actually exist::
951 sage: set_random_seed()
952 sage: n_max = max(2, JordanSpinEJA._max_random_instance_size())
953 sage: n = ZZ.random_element(2, n_max)
954 sage: J = JordanSpinEJA(n)
955 sage: y = J.random_element()
956 sage: while y == y.coefficient(0)*J.one():
957 ....: y = J.random_element()
958 sage: y0 = y.to_vector()[0]
959 sage: y_bar = y.to_vector()[1:]
960 sage: actual = y.minimal_polynomial()
961 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
962 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
963 sage: bool(actual == expected)
966 The minimal polynomial should always kill its element::
968 sage: set_random_seed()
969 sage: x = random_eja().random_element()
970 sage: p = x.minimal_polynomial()
971 sage: x.apply_univariate_polynomial(p)
974 The minimal polynomial is invariant under a change of basis,
975 and in particular, a re-scaling of the basis::
977 sage: set_random_seed()
978 sage: n_max = RealSymmetricEJA._max_random_instance_size()
979 sage: n = ZZ.random_element(1, n_max)
980 sage: J1 = RealSymmetricEJA(n)
981 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
982 sage: X = random_matrix(AA,n)
983 sage: X = X*X.transpose()
986 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
991 # We would generate a zero-dimensional subalgebra
992 # where the minimal polynomial would be constant.
993 # That might be correct, but only if *this* algebra
995 if not self
.parent().is_trivial():
996 # Pretty sure we know what the minimal polynomial of
997 # the zero operator is going to be. This ensures
998 # consistency of e.g. the polynomial variable returned
999 # in the "normal" case without us having to think about it.
1000 return self
.operator().minimal_polynomial()
1002 A
= self
.subalgebra_generated_by(orthonormalize_basis
=False)
1003 return A(self
).operator().minimal_polynomial()
1007 def to_matrix(self
):
1009 Return an (often more natural) representation of this element as a
1012 Every finite-dimensional Euclidean Jordan Algebra is a direct
1013 sum of five simple algebras, four of which comprise Hermitian
1014 matrices. This method returns a "natural" matrix
1015 representation of this element as either a Hermitian matrix or
1020 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1021 ....: QuaternionHermitianEJA)
1025 sage: J = ComplexHermitianEJA(3)
1028 sage: J.one().to_matrix()
1038 sage: J = QuaternionHermitianEJA(3)
1041 sage: J.one().to_matrix()
1042 [1 0 0 0 0 0 0 0 0 0 0 0]
1043 [0 1 0 0 0 0 0 0 0 0 0 0]
1044 [0 0 1 0 0 0 0 0 0 0 0 0]
1045 [0 0 0 1 0 0 0 0 0 0 0 0]
1046 [0 0 0 0 1 0 0 0 0 0 0 0]
1047 [0 0 0 0 0 1 0 0 0 0 0 0]
1048 [0 0 0 0 0 0 1 0 0 0 0 0]
1049 [0 0 0 0 0 0 0 1 0 0 0 0]
1050 [0 0 0 0 0 0 0 0 1 0 0 0]
1051 [0 0 0 0 0 0 0 0 0 1 0 0]
1052 [0 0 0 0 0 0 0 0 0 0 1 0]
1053 [0 0 0 0 0 0 0 0 0 0 0 1]
1055 B
= self
.parent().matrix_basis()
1056 W
= self
.parent().matrix_space()
1058 # This is just a manual "from_vector()", but of course
1059 # matrix spaces aren't vector spaces in sage, so they
1060 # don't have a from_vector() method.
1061 return W
.linear_combination( zip(B
, self
.to_vector()) )
1066 The norm of this element with respect to :meth:`inner_product`.
1070 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1075 sage: J = HadamardEJA(2)
1076 sage: x = sum(J.gens())
1082 sage: J = JordanSpinEJA(4)
1083 sage: x = sum(J.gens())
1088 return self
.inner_product(self
).sqrt()
1093 Return the left-multiplication-by-this-element
1094 operator on the ambient algebra.
1098 sage: from mjo.eja.eja_algebra import random_eja
1102 sage: set_random_seed()
1103 sage: J = random_eja()
1104 sage: x,y = J.random_elements(2)
1105 sage: x.operator()(y) == x*y
1107 sage: y.operator()(x) == x*y
1112 left_mult_by_self
= lambda y
: self
*y
1113 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1114 return FiniteDimensionalEuclideanJordanAlgebraOperator(
1120 def quadratic_representation(self
, other
=None):
1122 Return the quadratic representation of this element.
1126 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1131 The explicit form in the spin factor algebra is given by
1132 Alizadeh's Example 11.12::
1134 sage: set_random_seed()
1135 sage: x = JordanSpinEJA.random_instance().random_element()
1136 sage: x_vec = x.to_vector()
1137 sage: Q = matrix.identity(x.base_ring(), 0)
1138 sage: n = x_vec.degree()
1141 ....: x_bar = x_vec[1:]
1142 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1143 ....: B = 2*x0*x_bar.row()
1144 ....: C = 2*x0*x_bar.column()
1145 ....: D = matrix.identity(x.base_ring(), n-1)
1146 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1147 ....: D = D + 2*x_bar.tensor_product(x_bar)
1148 ....: Q = matrix.block(2,2,[A,B,C,D])
1149 sage: Q == x.quadratic_representation().matrix()
1152 Test all of the properties from Theorem 11.2 in Alizadeh::
1154 sage: set_random_seed()
1155 sage: J = random_eja()
1156 sage: x,y = J.random_elements(2)
1157 sage: Lx = x.operator()
1158 sage: Lxx = (x*x).operator()
1159 sage: Qx = x.quadratic_representation()
1160 sage: Qy = y.quadratic_representation()
1161 sage: Qxy = x.quadratic_representation(y)
1162 sage: Qex = J.one().quadratic_representation(x)
1163 sage: n = ZZ.random_element(10)
1164 sage: Qxn = (x^n).quadratic_representation()
1168 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1171 Property 2 (multiply on the right for :trac:`28272`):
1173 sage: alpha = J.base_ring().random_element()
1174 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1179 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1182 sage: not x.is_invertible() or (
1185 ....: x.inverse().quadratic_representation() )
1188 sage: Qxy(J.one()) == x*y
1193 sage: not x.is_invertible() or (
1194 ....: x.quadratic_representation(x.inverse())*Qx
1195 ....: == Qx*x.quadratic_representation(x.inverse()) )
1198 sage: not x.is_invertible() or (
1199 ....: x.quadratic_representation(x.inverse())*Qx
1201 ....: 2*Lx*Qex - Qx )
1204 sage: 2*Lx*Qex - Qx == Lxx
1209 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1219 sage: not x.is_invertible() or (
1220 ....: Qx*x.inverse().operator() == Lx )
1225 sage: not x.operator_commutes_with(y) or (
1226 ....: Qx(y)^n == Qxn(y^n) )
1232 elif not other
in self
.parent():
1233 raise TypeError("'other' must live in the same algebra")
1236 M
= other
.operator()
1237 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1241 def spectral_decomposition(self
):
1243 Return the unique spectral decomposition of this element.
1247 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1248 element's left-multiplication-by operator to the subalgebra it
1249 generates. We then compute the spectral decomposition of that
1250 operator, and the spectral projectors we get back must be the
1251 left-multiplication-by operators for the idempotents we
1252 seek. Thus applying them to the identity element gives us those
1255 Since the eigenvalues are required to be distinct, we take
1256 the spectral decomposition of the zero element to be zero
1257 times the identity element of the algebra (which is idempotent,
1262 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1266 The spectral decomposition of the identity is ``1`` times itself,
1267 and the spectral decomposition of zero is ``0`` times the identity::
1269 sage: J = RealSymmetricEJA(3)
1272 sage: J.one().spectral_decomposition()
1274 sage: J.zero().spectral_decomposition()
1279 sage: J = RealSymmetricEJA(4)
1280 sage: x = sum(J.gens())
1281 sage: sd = x.spectral_decomposition()
1286 sage: c0.inner_product(c1) == 0
1288 sage: c0.is_idempotent()
1290 sage: c1.is_idempotent()
1292 sage: c0 + c1 == J.one()
1294 sage: l0*c0 + l1*c1 == x
1297 The spectral decomposition should work in subalgebras, too::
1299 sage: J = RealSymmetricEJA(4)
1300 sage: (e0, e1, e2, e3, e4, e5, e6, e7, e8, e9) = J.gens()
1301 sage: A = 2*e5 - 2*e8
1302 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1303 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1304 sage: (f0, f1, f2) = J1.gens()
1305 sage: f0.spectral_decomposition()
1309 A
= self
.subalgebra_generated_by(orthonormalize_basis
=True)
1311 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1312 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1315 def subalgebra_generated_by(self
, orthonormalize_basis
=False):
1317 Return the associative subalgebra of the parent EJA generated
1320 Since our parent algebra is unital, we want "subalgebra" to mean
1321 "unital subalgebra" as well; thus the subalgebra that an element
1322 generates will itself be a Euclidean Jordan algebra after
1323 restricting the algebra operations appropriately. This is the
1324 subalgebra that Faraut and Korányi work with in section II.2, for
1329 sage: from mjo.eja.eja_algebra import random_eja
1333 This subalgebra, being composed of only powers, is associative::
1335 sage: set_random_seed()
1336 sage: x0 = random_eja().random_element()
1337 sage: A = x0.subalgebra_generated_by()
1338 sage: x,y,z = A.random_elements(3)
1339 sage: (x*y)*z == x*(y*z)
1342 Squaring in the subalgebra should work the same as in
1345 sage: set_random_seed()
1346 sage: x = random_eja().random_element()
1347 sage: A = x.subalgebra_generated_by()
1348 sage: A(x^2) == A(x)*A(x)
1351 By definition, the subalgebra generated by the zero element is
1352 the one-dimensional algebra generated by the identity
1353 element... unless the original algebra was trivial, in which
1354 case the subalgebra is trivial too::
1356 sage: set_random_seed()
1357 sage: A = random_eja().zero().subalgebra_generated_by()
1358 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1362 from mjo
.eja
.eja_element_subalgebra
import FiniteDimensionalEuclideanJordanElementSubalgebra
1363 return FiniteDimensionalEuclideanJordanElementSubalgebra(self
, orthonormalize_basis
)
1366 def subalgebra_idempotent(self
):
1368 Find an idempotent in the associative subalgebra I generate
1369 using Proposition 2.3.5 in Baes.
1373 sage: from mjo.eja.eja_algebra import random_eja
1377 Ensure that we can find an idempotent in a non-trivial algebra
1378 where there are non-nilpotent elements, or that we get the dumb
1379 solution in the trivial algebra::
1381 sage: set_random_seed()
1382 sage: J = random_eja()
1383 sage: x = J.random_element()
1384 sage: while x.is_nilpotent() and not J.is_trivial():
1385 ....: x = J.random_element()
1386 sage: c = x.subalgebra_idempotent()
1391 if self
.parent().is_trivial():
1394 if self
.is_nilpotent():
1395 raise ValueError("this only works with non-nilpotent elements!")
1397 J
= self
.subalgebra_generated_by()
1400 # The image of the matrix of left-u^m-multiplication
1401 # will be minimal for some natural number s...
1403 minimal_dim
= J
.dimension()
1404 for i
in range(1, minimal_dim
):
1405 this_dim
= (u
**i
).operator().matrix().image().dimension()
1406 if this_dim
< minimal_dim
:
1407 minimal_dim
= this_dim
1410 # Now minimal_matrix should correspond to the smallest
1411 # non-zero subspace in Baes's (or really, Koecher's)
1414 # However, we need to restrict the matrix to work on the
1415 # subspace... or do we? Can't we just solve, knowing that
1416 # A(c) = u^(s+1) should have a solution in the big space,
1419 # Beware, solve_right() means that we're using COLUMN vectors.
1420 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1422 A
= u_next
.operator().matrix()
1423 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1425 # Now c is the idempotent we want, but it still lives in the subalgebra.
1426 return c
.superalgebra_element()
1431 Return my trace, the sum of my eigenvalues.
1433 In a trivial algebra, however you want to look at it, the trace is
1434 an empty sum for which we declare the result to be zero.
1438 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1445 sage: J = TrivialEJA()
1446 sage: J.zero().trace()
1450 sage: J = JordanSpinEJA(3)
1451 sage: x = sum(J.gens())
1457 sage: J = HadamardEJA(5)
1458 sage: J.one().trace()
1463 The trace of an element is a real number::
1465 sage: set_random_seed()
1466 sage: J = random_eja()
1467 sage: J.random_element().trace() in RLF
1475 # Special case for the trivial algebra where
1476 # the trace is an empty sum.
1477 return P
.base_ring().zero()
1479 p
= P
._charpoly
_coefficients
()[r
-1]
1480 # The _charpoly_coeff function already adds the factor of
1481 # -1 to ensure that _charpoly_coeff(r-1) is really what
1482 # appears in front of t^{r-1} in the charpoly. However,
1483 # we want the negative of THAT for the trace.
1484 return -p(*self
.to_vector())
1487 def trace_inner_product(self
, other
):
1489 Return the trace inner product of myself and ``other``.
1493 sage: from mjo.eja.eja_algebra import random_eja
1497 The trace inner product is commutative, bilinear, and associative::
1499 sage: set_random_seed()
1500 sage: J = random_eja()
1501 sage: x,y,z = J.random_elements(3)
1503 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1506 sage: a = J.base_ring().random_element();
1507 sage: actual = (a*(x+z)).trace_inner_product(y)
1508 sage: expected = ( a*x.trace_inner_product(y) +
1509 ....: a*z.trace_inner_product(y) )
1510 sage: actual == expected
1512 sage: actual = x.trace_inner_product(a*(y+z))
1513 sage: expected = ( a*x.trace_inner_product(y) +
1514 ....: a*x.trace_inner_product(z) )
1515 sage: actual == expected
1518 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1522 if not other
in self
.parent():
1523 raise TypeError("'other' must live in the same algebra")
1525 return (self
*other
).trace()
1528 def trace_norm(self
):
1530 The norm of this element with respect to :meth:`trace_inner_product`.
1534 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1539 sage: J = HadamardEJA(2)
1540 sage: x = sum(J.gens())
1541 sage: x.trace_norm()
1546 sage: J = JordanSpinEJA(4)
1547 sage: x = sum(J.gens())
1548 sage: x.trace_norm()
1552 return self
.trace_inner_product(self
).sqrt()