1 from sage
.matrix
.constructor
import matrix
2 from sage
.misc
.cachefunc
import cached_method
3 from sage
.modules
.free_module
import VectorSpace
4 from sage
.modules
.with_basis
.indexed_element
import IndexedFreeModuleElement
6 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
7 from mjo
.eja
.eja_utils
import _mat2vec
9 class FiniteDimensionalEJAElement(IndexedFreeModuleElement
):
11 An element of a Euclidean Jordan algebra.
16 Oh man, I should not be doing this. This hides the "disabled"
17 methods ``left_matrix`` and ``matrix`` from introspection;
18 in particular it removes them from tab-completion.
20 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
28 Return ``self`` raised to the power ``n``.
30 Jordan algebras are always power-associative; see for
31 example Faraut and Korányi, Proposition II.1.2 (ii).
33 We have to override this because our superclass uses row
34 vectors instead of column vectors! We, on the other hand,
35 assume column vectors everywhere.
39 sage: from mjo.eja.eja_algebra import random_eja
43 The definition of `x^2` is the unambiguous `x*x`::
45 sage: set_random_seed()
46 sage: x = random_eja().random_element()
50 A few examples of power-associativity::
52 sage: set_random_seed()
53 sage: x = random_eja().random_element()
54 sage: x*(x*x)*(x*x) == x^5
56 sage: (x*x)*(x*x*x) == x^5
59 We also know that powers operator-commute (Koecher, Chapter
62 sage: set_random_seed()
63 sage: x = random_eja().random_element()
64 sage: m = ZZ.random_element(0,10)
65 sage: n = ZZ.random_element(0,10)
66 sage: Lxm = (x^m).operator()
67 sage: Lxn = (x^n).operator()
68 sage: Lxm*Lxn == Lxn*Lxm
73 return self
.parent().one()
77 return (self
**(n
-1))*self
80 def apply_univariate_polynomial(self
, p
):
82 Apply the univariate polynomial ``p`` to this element.
84 A priori, SageMath won't allow us to apply a univariate
85 polynomial to an element of an EJA, because we don't know
86 that EJAs are rings (they are usually not associative). Of
87 course, we know that EJAs are power-associative, so the
88 operation is ultimately kosher. This function sidesteps
89 the CAS to get the answer we want and expect.
93 sage: from mjo.eja.eja_algebra import (HadamardEJA,
98 sage: R = PolynomialRing(QQ, 't')
100 sage: p = t^4 - t^3 + 5*t - 2
101 sage: J = HadamardEJA(5)
102 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
107 We should always get back an element of the algebra::
109 sage: set_random_seed()
110 sage: p = PolynomialRing(AA, 't').random_element()
111 sage: J = random_eja()
112 sage: x = J.random_element()
113 sage: x.apply_univariate_polynomial(p) in J
117 if len(p
.variables()) > 1:
118 raise ValueError("not a univariate polynomial")
121 # Convert the coeficcients to the parent's base ring,
122 # because a priori they might live in an (unnecessarily)
123 # larger ring for which P.sum() would fail below.
124 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
125 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
128 def characteristic_polynomial(self
):
130 Return the characteristic polynomial of this element.
134 sage: from mjo.eja.eja_algebra import HadamardEJA
138 The rank of `R^3` is three, and the minimal polynomial of
139 the identity element is `(t-1)` from which it follows that
140 the characteristic polynomial should be `(t-1)^3`::
142 sage: J = HadamardEJA(3)
143 sage: J.one().characteristic_polynomial()
144 t^3 - 3*t^2 + 3*t - 1
146 Likewise, the characteristic of the zero element in the
147 rank-three algebra `R^{n}` should be `t^{3}`::
149 sage: J = HadamardEJA(3)
150 sage: J.zero().characteristic_polynomial()
155 The characteristic polynomial of an element should evaluate
156 to zero on that element::
158 sage: set_random_seed()
159 sage: x = HadamardEJA(3).random_element()
160 sage: p = x.characteristic_polynomial()
161 sage: x.apply_univariate_polynomial(p)
164 The characteristic polynomials of the zero and unit elements
165 should be what we think they are in a subalgebra, too::
167 sage: J = HadamardEJA(3)
168 sage: p1 = J.one().characteristic_polynomial()
169 sage: q1 = J.zero().characteristic_polynomial()
170 sage: e0,e1,e2 = J.gens()
171 sage: A = (e0 + 2*e1 + 3*e2).subalgebra_generated_by() # dim 3
172 sage: p2 = A.one().characteristic_polynomial()
173 sage: q2 = A.zero().characteristic_polynomial()
180 p
= self
.parent().characteristic_polynomial_of()
181 return p(*self
.to_vector())
184 def inner_product(self
, other
):
186 Return the parent algebra's inner product of myself and ``other``.
190 sage: from mjo.eja.eja_algebra import (
191 ....: ComplexHermitianEJA,
193 ....: QuaternionHermitianEJA,
194 ....: RealSymmetricEJA,
199 The inner product in the Jordan spin algebra is the usual
200 inner product on `R^n` (this example only works because the
201 basis for the Jordan algebra is the standard basis in `R^n`)::
203 sage: J = JordanSpinEJA(3)
204 sage: x = vector(QQ,[1,2,3])
205 sage: y = vector(QQ,[4,5,6])
206 sage: x.inner_product(y)
208 sage: J.from_vector(x).inner_product(J.from_vector(y))
211 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
212 multiplication is the usual matrix multiplication in `S^n`,
213 so the inner product of the identity matrix with itself
216 sage: J = RealSymmetricEJA(3)
217 sage: J.one().inner_product(J.one())
220 Likewise, the inner product on `C^n` is `<X,Y> =
221 Re(trace(X*Y))`, where we must necessarily take the real
222 part because the product of Hermitian matrices may not be
225 sage: J = ComplexHermitianEJA(3)
226 sage: J.one().inner_product(J.one())
229 Ditto for the quaternions::
231 sage: J = QuaternionHermitianEJA(2)
232 sage: J.one().inner_product(J.one())
237 Ensure that we can always compute an inner product, and that
238 it gives us back a real number::
240 sage: set_random_seed()
241 sage: J = random_eja()
242 sage: x,y = J.random_elements(2)
243 sage: x.inner_product(y) in RLF
249 raise TypeError("'other' must live in the same algebra")
251 return P
.inner_product(self
, other
)
254 def operator_commutes_with(self
, other
):
256 Return whether or not this element operator-commutes
261 sage: from mjo.eja.eja_algebra import random_eja
265 The definition of a Jordan algebra says that any element
266 operator-commutes with its square::
268 sage: set_random_seed()
269 sage: x = random_eja().random_element()
270 sage: x.operator_commutes_with(x^2)
275 Test Lemma 1 from Chapter III of Koecher::
277 sage: set_random_seed()
278 sage: u,v = random_eja().random_elements(2)
279 sage: lhs = u.operator_commutes_with(u*v)
280 sage: rhs = v.operator_commutes_with(u^2)
284 Test the first polarization identity from my notes, Koecher
285 Chapter III, or from Baes (2.3)::
287 sage: set_random_seed()
288 sage: x,y = random_eja().random_elements(2)
289 sage: Lx = x.operator()
290 sage: Ly = y.operator()
291 sage: Lxx = (x*x).operator()
292 sage: Lxy = (x*y).operator()
293 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
296 Test the second polarization identity from my notes or from
299 sage: set_random_seed()
300 sage: x,y,z = random_eja().random_elements(3)
301 sage: Lx = x.operator()
302 sage: Ly = y.operator()
303 sage: Lz = z.operator()
304 sage: Lzy = (z*y).operator()
305 sage: Lxy = (x*y).operator()
306 sage: Lxz = (x*z).operator()
307 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
310 Test the third polarization identity from my notes or from
313 sage: set_random_seed()
314 sage: u,y,z = random_eja().random_elements(3)
315 sage: Lu = u.operator()
316 sage: Ly = y.operator()
317 sage: Lz = z.operator()
318 sage: Lzy = (z*y).operator()
319 sage: Luy = (u*y).operator()
320 sage: Luz = (u*z).operator()
321 sage: Luyz = (u*(y*z)).operator()
322 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
323 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
324 sage: bool(lhs == rhs)
328 if not other
in self
.parent():
329 raise TypeError("'other' must live in the same algebra")
338 Return my determinant, the product of my eigenvalues.
342 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
344 ....: RealSymmetricEJA,
345 ....: ComplexHermitianEJA,
350 sage: J = JordanSpinEJA(2)
351 sage: e0,e1 = J.gens()
352 sage: x = sum( J.gens() )
358 sage: J = JordanSpinEJA(3)
359 sage: e0,e1,e2 = J.gens()
360 sage: x = sum( J.gens() )
364 The determinant of the sole element in the rank-zero trivial
365 algebra is ``1``, by three paths of reasoning. First, its
366 characteristic polynomial is a constant ``1``, so the constant
367 term in that polynomial is ``1``. Second, the characteristic
368 polynomial evaluated at zero is again ``1``. And finally, the
369 (empty) product of its eigenvalues is likewise just unity::
371 sage: J = TrivialEJA()
377 An element is invertible if and only if its determinant is
380 sage: set_random_seed()
381 sage: x = random_eja().random_element()
382 sage: x.is_invertible() == (x.det() != 0)
385 Ensure that the determinant is multiplicative on an associative
386 subalgebra as in Faraut and Korányi's Proposition II.2.2::
388 sage: set_random_seed()
389 sage: J = random_eja().random_element().subalgebra_generated_by()
390 sage: x,y = J.random_elements(2)
391 sage: (x*y).det() == x.det()*y.det()
394 The determinant in matrix algebras is just the usual determinant::
396 sage: set_random_seed()
397 sage: X = matrix.random(QQ,3)
399 sage: J1 = RealSymmetricEJA(3)
400 sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
401 sage: expected = X.det()
402 sage: actual1 = J1(X).det()
403 sage: actual2 = J2(X).det()
404 sage: actual1 == expected
406 sage: actual2 == expected
411 sage: set_random_seed()
412 sage: J1 = ComplexHermitianEJA(2)
413 sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
414 sage: X = matrix.random(GaussianIntegers(), 2)
416 sage: expected = AA(X.det())
417 sage: actual1 = J1(J1.real_embed(X)).det()
418 sage: actual2 = J2(J2.real_embed(X)).det()
419 sage: expected == actual1
421 sage: expected == actual2
429 # Special case, since we don't get the a0=1
430 # coefficient when the rank of the algebra
432 return P
.base_ring().one()
434 p
= P
._charpoly
_coefficients
()[0]
435 # The _charpoly_coeff function already adds the factor of -1
436 # to ensure that _charpoly_coefficients()[0] is really what
437 # appears in front of t^{0} in the charpoly. However, we want
438 # (-1)^r times THAT for the determinant.
439 return ((-1)**r
)*p(*self
.to_vector())
445 Return the Jordan-multiplicative inverse of this element.
449 In general we appeal to the quadratic representation as in
450 Koecher's Theorem 12 in Chapter III, Section 5. But if the
451 parent algebra's "characteristic polynomial of" coefficients
452 happen to be cached, then we use Proposition II.2.4 in Faraut
453 and Korányi which gives a formula for the inverse based on the
454 characteristic polynomial and the Cayley-Hamilton theorem for
455 Euclidean Jordan algebras::
459 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
465 The inverse in the spin factor algebra is given in Alizadeh's
468 sage: set_random_seed()
469 sage: J = JordanSpinEJA.random_instance()
470 sage: x = J.random_element()
471 sage: while not x.is_invertible():
472 ....: x = J.random_element()
473 sage: x_vec = x.to_vector()
475 sage: x_bar = x_vec[1:]
476 sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
477 sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
478 sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
479 sage: x.inverse() == J.from_vector(x_inverse)
482 Trying to invert a non-invertible element throws an error:
484 sage: JordanSpinEJA(3).zero().inverse()
485 Traceback (most recent call last):
487 ZeroDivisionError: element is not invertible
491 The identity element is its own inverse::
493 sage: set_random_seed()
494 sage: J = random_eja()
495 sage: J.one().inverse() == J.one()
498 If an element has an inverse, it acts like one::
500 sage: set_random_seed()
501 sage: J = random_eja()
502 sage: x = J.random_element()
503 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
506 The inverse of the inverse is what we started with::
508 sage: set_random_seed()
509 sage: J = random_eja()
510 sage: x = J.random_element()
511 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
514 Proposition II.2.3 in Faraut and Korányi says that the inverse
515 of an element is the inverse of its left-multiplication operator
516 applied to the algebra's identity, when that inverse exists::
518 sage: set_random_seed()
519 sage: J = random_eja()
520 sage: x = J.random_element()
521 sage: (not x.operator().is_invertible()) or (
522 ....: x.operator().inverse()(J.one()) == x.inverse() )
525 Check that the fast (cached) and slow algorithms give the same
528 sage: set_random_seed() # long time
529 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
530 sage: x = J.random_element() # long time
531 sage: while not x.is_invertible(): # long time
532 ....: x = J.random_element() # long time
533 sage: slow = x.inverse() # long time
534 sage: _ = J._charpoly_coefficients() # long time
535 sage: fast = x.inverse() # long time
536 sage: slow == fast # long time
539 not_invertible_msg
= "element is not invertible"
540 if self
.parent()._charpoly
_coefficients
.is_in_cache():
541 # We can invert using our charpoly if it will be fast to
542 # compute. If the coefficients are cached, our rank had
544 if self
.det().is_zero():
545 raise ZeroDivisionError(not_invertible_msg
)
546 r
= self
.parent().rank()
547 a
= self
.characteristic_polynomial().coefficients(sparse
=False)
548 return (-1)**(r
+1)*sum(a
[i
+1]*self
**i
for i
in range(r
))/self
.det()
551 inv
= (~self
.quadratic_representation())(self
)
552 self
.is_invertible
.set_cache(True)
554 except ZeroDivisionError:
555 self
.is_invertible
.set_cache(False)
556 raise ZeroDivisionError(not_invertible_msg
)
560 def is_invertible(self
):
562 Return whether or not this element is invertible.
566 If computing my determinant will be fast, we do so and compare
567 with zero (Proposition II.2.4 in Faraut and
568 Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
569 reduces the problem to the invertibility of my quadratic
574 sage: from mjo.eja.eja_algebra import random_eja
578 The identity element is always invertible::
580 sage: set_random_seed()
581 sage: J = random_eja()
582 sage: J.one().is_invertible()
585 The zero element is never invertible in a non-trivial algebra::
587 sage: set_random_seed()
588 sage: J = random_eja()
589 sage: (not J.is_trivial()) and J.zero().is_invertible()
592 Test that the fast (cached) and slow algorithms give the same
595 sage: set_random_seed() # long time
596 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
597 sage: x = J.random_element() # long time
598 sage: slow = x.is_invertible() # long time
599 sage: _ = J._charpoly_coefficients() # long time
600 sage: fast = x.is_invertible() # long time
601 sage: slow == fast # long time
605 if self
.parent().is_trivial():
610 if self
.parent()._charpoly
_coefficients
.is_in_cache():
611 # The determinant will be quicker than inverting the
612 # quadratic representation, most likely.
613 return (not self
.det().is_zero())
615 # The easiest way to determine if I'm invertible is to try.
617 inv
= (~self
.quadratic_representation())(self
)
618 self
.inverse
.set_cache(inv
)
620 except ZeroDivisionError:
624 def is_primitive_idempotent(self
):
626 Return whether or not this element is a primitive (or minimal)
629 A primitive idempotent is a non-zero idempotent that is not
630 the sum of two other non-zero idempotents. Remark 2.7.15 in
631 Baes shows that this is what he refers to as a "minimal
634 An element of a Euclidean Jordan algebra is a minimal idempotent
635 if it :meth:`is_idempotent` and if its Peirce subalgebra
636 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
641 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
642 ....: RealSymmetricEJA,
648 This method is sloooooow.
652 The spectral decomposition of a non-regular element should always
653 contain at least one non-minimal idempotent::
655 sage: J = RealSymmetricEJA(3)
656 sage: x = sum(J.gens())
659 sage: [ c.is_primitive_idempotent()
660 ....: for (l,c) in x.spectral_decomposition() ]
663 On the other hand, the spectral decomposition of a regular
664 element should always be in terms of minimal idempotents::
666 sage: J = JordanSpinEJA(4)
667 sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) )
670 sage: [ c.is_primitive_idempotent()
671 ....: for (l,c) in x.spectral_decomposition() ]
676 The identity element is minimal only in an EJA of rank one::
678 sage: set_random_seed()
679 sage: J = random_eja()
680 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
683 A non-idempotent cannot be a minimal idempotent::
685 sage: set_random_seed()
686 sage: J = JordanSpinEJA(4)
687 sage: x = J.random_element()
688 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
691 Proposition 2.7.19 in Baes says that an element is a minimal
692 idempotent if and only if it's idempotent with trace equal to
695 sage: set_random_seed()
696 sage: J = JordanSpinEJA(4)
697 sage: x = J.random_element()
698 sage: expected = (x.is_idempotent() and x.trace() == 1)
699 sage: actual = x.is_primitive_idempotent()
700 sage: actual == expected
703 Primitive idempotents must be non-zero::
705 sage: set_random_seed()
706 sage: J = random_eja()
707 sage: J.zero().is_idempotent()
709 sage: J.zero().is_primitive_idempotent()
712 As a consequence of the fact that primitive idempotents must
713 be non-zero, there are no primitive idempotents in a trivial
714 Euclidean Jordan algebra::
716 sage: J = TrivialEJA()
717 sage: J.one().is_idempotent()
719 sage: J.one().is_primitive_idempotent()
723 if not self
.is_idempotent():
729 (_
,_
,J1
) = self
.parent().peirce_decomposition(self
)
730 return (J1
.dimension() == 1)
733 def is_nilpotent(self
):
735 Return whether or not some power of this element is zero.
739 We use Theorem 5 in Chapter III of Koecher, which says that
740 an element ``x`` is nilpotent if and only if ``x.operator()``
741 is nilpotent. And it is a basic fact of linear algebra that
742 an operator on an `n`-dimensional space is nilpotent if and
743 only if, when raised to the `n`th power, it equals the zero
744 operator (for example, see Axler Corollary 8.8).
748 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
753 sage: J = JordanSpinEJA(3)
754 sage: x = sum(J.gens())
755 sage: x.is_nilpotent()
760 The identity element is never nilpotent, except in a trivial EJA::
762 sage: set_random_seed()
763 sage: J = random_eja()
764 sage: J.one().is_nilpotent() and not J.is_trivial()
767 The additive identity is always nilpotent::
769 sage: set_random_seed()
770 sage: random_eja().zero().is_nilpotent()
775 zero_operator
= P
.zero().operator()
776 return self
.operator()**P
.dimension() == zero_operator
779 def is_regular(self
):
781 Return whether or not this is a regular element.
785 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
790 The identity element always has degree one, but any element
791 linearly-independent from it is regular::
793 sage: J = JordanSpinEJA(5)
794 sage: J.one().is_regular()
796 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
797 sage: for x in J.gens():
798 ....: (J.one() + x).is_regular()
807 The zero element should never be regular, unless the parent
808 algebra has dimension less than or equal to one::
810 sage: set_random_seed()
811 sage: J = random_eja()
812 sage: J.dimension() <= 1 or not J.zero().is_regular()
815 The unit element isn't regular unless the algebra happens to
816 consist of only its scalar multiples::
818 sage: set_random_seed()
819 sage: J = random_eja()
820 sage: J.dimension() <= 1 or not J.one().is_regular()
824 return self
.degree() == self
.parent().rank()
829 Return the degree of this element, which is defined to be
830 the degree of its minimal polynomial.
838 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
843 sage: J = JordanSpinEJA(4)
844 sage: J.one().degree()
846 sage: e0,e1,e2,e3 = J.gens()
847 sage: (e0 - e1).degree()
850 In the spin factor algebra (of rank two), all elements that
851 aren't multiples of the identity are regular::
853 sage: set_random_seed()
854 sage: J = JordanSpinEJA.random_instance()
855 sage: n = J.dimension()
856 sage: x = J.random_element()
857 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
862 The zero and unit elements are both of degree one in nontrivial
865 sage: set_random_seed()
866 sage: J = random_eja()
867 sage: d = J.zero().degree()
868 sage: (J.is_trivial() and d == 0) or d == 1
870 sage: d = J.one().degree()
871 sage: (J.is_trivial() and d == 0) or d == 1
874 Our implementation agrees with the definition::
876 sage: set_random_seed()
877 sage: x = random_eja().random_element()
878 sage: x.degree() == x.minimal_polynomial().degree()
882 n
= self
.parent().dimension()
885 # The minimal polynomial is an empty product, i.e. the
886 # constant polynomial "1" having degree zero.
889 # The minimal polynomial of zero in a nontrivial algebra
890 # is "t", and is of degree one.
893 # If this is a nonzero element of a nontrivial algebra, it
894 # has degree at least one. It follows that, in an algebra
895 # of dimension one, the degree must be actually one.
898 # BEWARE: The subalgebra_generated_by() method uses the result
899 # of this method to construct a basis for the subalgebra. That
900 # means, in particular, that we cannot implement this method
901 # as ``self.subalgebra_generated_by().dimension()``.
903 # Algorithm: keep appending (vector representations of) powers
904 # self as rows to a matrix and echelonizing it. When its rank
905 # stops increasing, we've reached a redundancy.
907 # Given the special cases above, we can assume that "self" is
908 # nonzero, the algebra is nontrivial, and that its dimension
910 M
= matrix([(self
.parent().one()).to_vector()])
913 # Specifying the row-reduction algorithm can e.g. help over
914 # AA because it avoids the RecursionError that gets thrown
915 # when we have to look too hard for a root.
917 # Beware: QQ supports an entirely different set of "algorithm"
918 # keywords than do AA and RR.
920 from sage
.rings
.all
import QQ
921 if self
.parent().base_ring() is not QQ
:
922 algo
= "scaled_partial_pivoting"
925 M
= matrix(M
.rows() + [(self
**d
).to_vector()])
928 if new_rank
== old_rank
:
937 def left_matrix(self
):
939 Our parent class defines ``left_matrix`` and ``matrix``
940 methods whose names are misleading. We don't want them.
942 raise NotImplementedError("use operator().matrix() instead")
947 def minimal_polynomial(self
):
949 Return the minimal polynomial of this element,
950 as a function of the variable `t`.
954 We restrict ourselves to the associative subalgebra
955 generated by this element, and then return the minimal
956 polynomial of this element's operator matrix (in that
957 subalgebra). This works by Baes Proposition 2.3.16.
961 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
962 ....: RealSymmetricEJA,
968 Keeping in mind that the polynomial ``1`` evaluates the identity
969 element (also the zero element) of the trivial algebra, it is clear
970 that the polynomial ``1`` is the minimal polynomial of the only
971 element in a trivial algebra::
973 sage: J = TrivialEJA()
974 sage: J.one().minimal_polynomial()
976 sage: J.zero().minimal_polynomial()
981 The minimal polynomial of the identity and zero elements are
982 always the same, except in trivial algebras where the minimal
983 polynomial of the unit/zero element is ``1``::
985 sage: set_random_seed()
986 sage: J = random_eja()
987 sage: mu = J.one().minimal_polynomial()
988 sage: t = mu.parent().gen()
989 sage: mu + int(J.is_trivial())*(t-2)
991 sage: mu = J.zero().minimal_polynomial()
992 sage: t = mu.parent().gen()
993 sage: mu + int(J.is_trivial())*(t-1)
996 The degree of an element is (by one definition) the degree
997 of its minimal polynomial::
999 sage: set_random_seed()
1000 sage: x = random_eja().random_element()
1001 sage: x.degree() == x.minimal_polynomial().degree()
1004 The minimal polynomial and the characteristic polynomial coincide
1005 and are known (see Alizadeh, Example 11.11) for all elements of
1006 the spin factor algebra that aren't scalar multiples of the
1007 identity. We require the dimension of the algebra to be at least
1008 two here so that said elements actually exist::
1010 sage: set_random_seed()
1011 sage: n_max = max(2, JordanSpinEJA._max_random_instance_size())
1012 sage: n = ZZ.random_element(2, n_max)
1013 sage: J = JordanSpinEJA(n)
1014 sage: y = J.random_element()
1015 sage: while y == y.coefficient(0)*J.one():
1016 ....: y = J.random_element()
1017 sage: y0 = y.to_vector()[0]
1018 sage: y_bar = y.to_vector()[1:]
1019 sage: actual = y.minimal_polynomial()
1020 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1021 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1022 sage: bool(actual == expected)
1025 The minimal polynomial should always kill its element::
1027 sage: set_random_seed()
1028 sage: x = random_eja().random_element()
1029 sage: p = x.minimal_polynomial()
1030 sage: x.apply_univariate_polynomial(p)
1033 The minimal polynomial is invariant under a change of basis,
1034 and in particular, a re-scaling of the basis::
1036 sage: set_random_seed()
1037 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1038 sage: n = ZZ.random_element(1, n_max)
1039 sage: J1 = RealSymmetricEJA(n)
1040 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
1041 sage: X = random_matrix(AA,n)
1042 sage: X = X*X.transpose()
1045 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
1050 # We would generate a zero-dimensional subalgebra
1051 # where the minimal polynomial would be constant.
1052 # That might be correct, but only if *this* algebra
1054 if not self
.parent().is_trivial():
1055 # Pretty sure we know what the minimal polynomial of
1056 # the zero operator is going to be. This ensures
1057 # consistency of e.g. the polynomial variable returned
1058 # in the "normal" case without us having to think about it.
1059 return self
.operator().minimal_polynomial()
1061 A
= self
.subalgebra_generated_by(orthonormalize
=False)
1062 return A(self
).operator().minimal_polynomial()
1066 def to_matrix(self
):
1068 Return an (often more natural) representation of this element as a
1071 Every finite-dimensional Euclidean Jordan Algebra is a direct
1072 sum of five simple algebras, four of which comprise Hermitian
1073 matrices. This method returns a "natural" matrix
1074 representation of this element as either a Hermitian matrix or
1079 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1080 ....: QuaternionHermitianEJA)
1084 sage: J = ComplexHermitianEJA(3)
1087 sage: J.one().to_matrix()
1097 sage: J = QuaternionHermitianEJA(2)
1100 sage: J.one().to_matrix()
1111 B
= self
.parent().matrix_basis()
1112 W
= self
.parent().matrix_space()
1114 # This is just a manual "from_vector()", but of course
1115 # matrix spaces aren't vector spaces in sage, so they
1116 # don't have a from_vector() method.
1117 return W
.linear_combination( zip(B
, self
.to_vector()) )
1122 The norm of this element with respect to :meth:`inner_product`.
1126 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1131 sage: J = HadamardEJA(2)
1132 sage: x = sum(J.gens())
1138 sage: J = JordanSpinEJA(4)
1139 sage: x = sum(J.gens())
1144 return self
.inner_product(self
).sqrt()
1149 Return the left-multiplication-by-this-element
1150 operator on the ambient algebra.
1154 sage: from mjo.eja.eja_algebra import random_eja
1158 sage: set_random_seed()
1159 sage: J = random_eja()
1160 sage: x,y = J.random_elements(2)
1161 sage: x.operator()(y) == x*y
1163 sage: y.operator()(x) == x*y
1168 left_mult_by_self
= lambda y
: self
*y
1169 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1170 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1173 def quadratic_representation(self
, other
=None):
1175 Return the quadratic representation of this element.
1179 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1184 The explicit form in the spin factor algebra is given by
1185 Alizadeh's Example 11.12::
1187 sage: set_random_seed()
1188 sage: x = JordanSpinEJA.random_instance().random_element()
1189 sage: x_vec = x.to_vector()
1190 sage: Q = matrix.identity(x.base_ring(), 0)
1191 sage: n = x_vec.degree()
1194 ....: x_bar = x_vec[1:]
1195 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1196 ....: B = 2*x0*x_bar.row()
1197 ....: C = 2*x0*x_bar.column()
1198 ....: D = matrix.identity(x.base_ring(), n-1)
1199 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1200 ....: D = D + 2*x_bar.tensor_product(x_bar)
1201 ....: Q = matrix.block(2,2,[A,B,C,D])
1202 sage: Q == x.quadratic_representation().matrix()
1205 Test all of the properties from Theorem 11.2 in Alizadeh::
1207 sage: set_random_seed()
1208 sage: J = random_eja()
1209 sage: x,y = J.random_elements(2)
1210 sage: Lx = x.operator()
1211 sage: Lxx = (x*x).operator()
1212 sage: Qx = x.quadratic_representation()
1213 sage: Qy = y.quadratic_representation()
1214 sage: Qxy = x.quadratic_representation(y)
1215 sage: Qex = J.one().quadratic_representation(x)
1216 sage: n = ZZ.random_element(10)
1217 sage: Qxn = (x^n).quadratic_representation()
1221 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1224 Property 2 (multiply on the right for :trac:`28272`):
1226 sage: alpha = J.base_ring().random_element()
1227 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1232 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1235 sage: not x.is_invertible() or (
1238 ....: x.inverse().quadratic_representation() )
1241 sage: Qxy(J.one()) == x*y
1246 sage: not x.is_invertible() or (
1247 ....: x.quadratic_representation(x.inverse())*Qx
1248 ....: == Qx*x.quadratic_representation(x.inverse()) )
1251 sage: not x.is_invertible() or (
1252 ....: x.quadratic_representation(x.inverse())*Qx
1254 ....: 2*Lx*Qex - Qx )
1257 sage: 2*Lx*Qex - Qx == Lxx
1262 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1272 sage: not x.is_invertible() or (
1273 ....: Qx*x.inverse().operator() == Lx )
1278 sage: not x.operator_commutes_with(y) or (
1279 ....: Qx(y)^n == Qxn(y^n) )
1285 elif not other
in self
.parent():
1286 raise TypeError("'other' must live in the same algebra")
1289 M
= other
.operator()
1290 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1294 def spectral_decomposition(self
):
1296 Return the unique spectral decomposition of this element.
1300 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1301 element's left-multiplication-by operator to the subalgebra it
1302 generates. We then compute the spectral decomposition of that
1303 operator, and the spectral projectors we get back must be the
1304 left-multiplication-by operators for the idempotents we
1305 seek. Thus applying them to the identity element gives us those
1308 Since the eigenvalues are required to be distinct, we take
1309 the spectral decomposition of the zero element to be zero
1310 times the identity element of the algebra (which is idempotent,
1315 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1319 The spectral decomposition of the identity is ``1`` times itself,
1320 and the spectral decomposition of zero is ``0`` times the identity::
1322 sage: J = RealSymmetricEJA(3)
1325 sage: J.one().spectral_decomposition()
1327 sage: J.zero().spectral_decomposition()
1332 sage: J = RealSymmetricEJA(4)
1333 sage: x = sum(J.gens())
1334 sage: sd = x.spectral_decomposition()
1339 sage: c0.inner_product(c1) == 0
1341 sage: c0.is_idempotent()
1343 sage: c1.is_idempotent()
1345 sage: c0 + c1 == J.one()
1347 sage: l0*c0 + l1*c1 == x
1350 The spectral decomposition should work in subalgebras, too::
1352 sage: J = RealSymmetricEJA(4)
1353 sage: (e0, e1, e2, e3, e4, e5, e6, e7, e8, e9) = J.gens()
1354 sage: A = 2*e5 - 2*e8
1355 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1356 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1357 sage: (f0, f1, f2) = J1.gens()
1358 sage: f0.spectral_decomposition()
1362 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1364 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1365 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1368 def subalgebra_generated_by(self
, **kwargs
):
1370 Return the associative subalgebra of the parent EJA generated
1373 Since our parent algebra is unital, we want "subalgebra" to mean
1374 "unital subalgebra" as well; thus the subalgebra that an element
1375 generates will itself be a Euclidean Jordan algebra after
1376 restricting the algebra operations appropriately. This is the
1377 subalgebra that Faraut and Korányi work with in section II.2, for
1382 sage: from mjo.eja.eja_algebra import random_eja
1386 This subalgebra, being composed of only powers, is associative::
1388 sage: set_random_seed()
1389 sage: x0 = random_eja().random_element()
1390 sage: A = x0.subalgebra_generated_by()
1391 sage: x,y,z = A.random_elements(3)
1392 sage: (x*y)*z == x*(y*z)
1395 Squaring in the subalgebra should work the same as in
1398 sage: set_random_seed()
1399 sage: x = random_eja().random_element()
1400 sage: A = x.subalgebra_generated_by()
1401 sage: A(x^2) == A(x)*A(x)
1404 By definition, the subalgebra generated by the zero element is
1405 the one-dimensional algebra generated by the identity
1406 element... unless the original algebra was trivial, in which
1407 case the subalgebra is trivial too::
1409 sage: set_random_seed()
1410 sage: A = random_eja().zero().subalgebra_generated_by()
1411 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1415 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1416 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1417 A
= FiniteDimensionalEJASubalgebra(self
.parent(),
1421 A
.one
.set_cache(A(self
.parent().one()))
1425 def subalgebra_idempotent(self
):
1427 Find an idempotent in the associative subalgebra I generate
1428 using Proposition 2.3.5 in Baes.
1432 sage: from mjo.eja.eja_algebra import random_eja
1436 Ensure that we can find an idempotent in a non-trivial algebra
1437 where there are non-nilpotent elements, or that we get the dumb
1438 solution in the trivial algebra::
1440 sage: set_random_seed()
1441 sage: J = random_eja()
1442 sage: x = J.random_element()
1443 sage: while x.is_nilpotent() and not J.is_trivial():
1444 ....: x = J.random_element()
1445 sage: c = x.subalgebra_idempotent()
1450 if self
.parent().is_trivial():
1453 if self
.is_nilpotent():
1454 raise ValueError("this only works with non-nilpotent elements!")
1456 J
= self
.subalgebra_generated_by()
1459 # The image of the matrix of left-u^m-multiplication
1460 # will be minimal for some natural number s...
1462 minimal_dim
= J
.dimension()
1463 for i
in range(1, minimal_dim
):
1464 this_dim
= (u
**i
).operator().matrix().image().dimension()
1465 if this_dim
< minimal_dim
:
1466 minimal_dim
= this_dim
1469 # Now minimal_matrix should correspond to the smallest
1470 # non-zero subspace in Baes's (or really, Koecher's)
1473 # However, we need to restrict the matrix to work on the
1474 # subspace... or do we? Can't we just solve, knowing that
1475 # A(c) = u^(s+1) should have a solution in the big space,
1478 # Beware, solve_right() means that we're using COLUMN vectors.
1479 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1481 A
= u_next
.operator().matrix()
1482 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1484 # Now c is the idempotent we want, but it still lives in the subalgebra.
1485 return c
.superalgebra_element()
1490 Return my trace, the sum of my eigenvalues.
1492 In a trivial algebra, however you want to look at it, the trace is
1493 an empty sum for which we declare the result to be zero.
1497 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1504 sage: J = TrivialEJA()
1505 sage: J.zero().trace()
1509 sage: J = JordanSpinEJA(3)
1510 sage: x = sum(J.gens())
1516 sage: J = HadamardEJA(5)
1517 sage: J.one().trace()
1522 The trace of an element is a real number::
1524 sage: set_random_seed()
1525 sage: J = random_eja()
1526 sage: J.random_element().trace() in RLF
1529 The trace is linear::
1531 sage: set_random_seed()
1532 sage: J = random_eja()
1533 sage: x,y = J.random_elements(2)
1534 sage: alpha = J.base_ring().random_element()
1535 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1543 # Special case for the trivial algebra where
1544 # the trace is an empty sum.
1545 return P
.base_ring().zero()
1547 p
= P
._charpoly
_coefficients
()[r
-1]
1548 # The _charpoly_coeff function already adds the factor of
1549 # -1 to ensure that _charpoly_coeff(r-1) is really what
1550 # appears in front of t^{r-1} in the charpoly. However,
1551 # we want the negative of THAT for the trace.
1552 return -p(*self
.to_vector())
1555 def trace_inner_product(self
, other
):
1557 Return the trace inner product of myself and ``other``.
1561 sage: from mjo.eja.eja_algebra import random_eja
1565 The trace inner product is commutative, bilinear, and associative::
1567 sage: set_random_seed()
1568 sage: J = random_eja()
1569 sage: x,y,z = J.random_elements(3)
1571 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1574 sage: a = J.base_ring().random_element();
1575 sage: actual = (a*(x+z)).trace_inner_product(y)
1576 sage: expected = ( a*x.trace_inner_product(y) +
1577 ....: a*z.trace_inner_product(y) )
1578 sage: actual == expected
1580 sage: actual = x.trace_inner_product(a*(y+z))
1581 sage: expected = ( a*x.trace_inner_product(y) +
1582 ....: a*x.trace_inner_product(z) )
1583 sage: actual == expected
1586 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1590 if not other
in self
.parent():
1591 raise TypeError("'other' must live in the same algebra")
1593 return (self
*other
).trace()
1596 def trace_norm(self
):
1598 The norm of this element with respect to :meth:`trace_inner_product`.
1602 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1607 sage: J = HadamardEJA(2)
1608 sage: x = sum(J.gens())
1609 sage: x.trace_norm()
1614 sage: J = JordanSpinEJA(4)
1615 sage: x = sum(J.gens())
1616 sage: x.trace_norm()
1620 return self
.trace_inner_product(self
).sqrt()