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()
1110 This also works in Cartesian product algebras::
1112 sage: J1 = HadamardEJA(1)
1113 sage: J2 = RealSymmetricEJA(2)
1114 sage: J = cartesian_product([J1,J2])
1115 sage: x = sum(J.gens())
1116 sage: x.to_matrix()[0]
1118 sage: x.to_matrix()[1]
1119 [ 1 0.7071067811865475?]
1120 [0.7071067811865475? 1]
1123 B
= self
.parent().matrix_basis()
1124 W
= self
.parent().matrix_space()
1126 if self
.parent()._matrix
_basis
_is
_cartesian
:
1127 # Aaaaand linear combinations don't work in Cartesian
1128 # product spaces, even though they provide a method
1130 pairs
= zip(B
, self
.to_vector())
1131 return sum( ( W(tuple(alpha
*b_i
for b_i
in b
))
1132 for (b
,alpha
) in pairs
),
1135 # This is just a manual "from_vector()", but of course
1136 # matrix spaces aren't vector spaces in sage, so they
1137 # don't have a from_vector() method.
1138 return W
.linear_combination( zip(B
, self
.to_vector()) )
1144 The norm of this element with respect to :meth:`inner_product`.
1148 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1153 sage: J = HadamardEJA(2)
1154 sage: x = sum(J.gens())
1160 sage: J = JordanSpinEJA(4)
1161 sage: x = sum(J.gens())
1166 return self
.inner_product(self
).sqrt()
1171 Return the left-multiplication-by-this-element
1172 operator on the ambient algebra.
1176 sage: from mjo.eja.eja_algebra import random_eja
1180 sage: set_random_seed()
1181 sage: J = random_eja()
1182 sage: x,y = J.random_elements(2)
1183 sage: x.operator()(y) == x*y
1185 sage: y.operator()(x) == x*y
1190 left_mult_by_self
= lambda y
: self
*y
1191 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1192 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1195 def quadratic_representation(self
, other
=None):
1197 Return the quadratic representation of this element.
1201 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1206 The explicit form in the spin factor algebra is given by
1207 Alizadeh's Example 11.12::
1209 sage: set_random_seed()
1210 sage: x = JordanSpinEJA.random_instance().random_element()
1211 sage: x_vec = x.to_vector()
1212 sage: Q = matrix.identity(x.base_ring(), 0)
1213 sage: n = x_vec.degree()
1216 ....: x_bar = x_vec[1:]
1217 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1218 ....: B = 2*x0*x_bar.row()
1219 ....: C = 2*x0*x_bar.column()
1220 ....: D = matrix.identity(x.base_ring(), n-1)
1221 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1222 ....: D = D + 2*x_bar.tensor_product(x_bar)
1223 ....: Q = matrix.block(2,2,[A,B,C,D])
1224 sage: Q == x.quadratic_representation().matrix()
1227 Test all of the properties from Theorem 11.2 in Alizadeh::
1229 sage: set_random_seed()
1230 sage: J = random_eja()
1231 sage: x,y = J.random_elements(2)
1232 sage: Lx = x.operator()
1233 sage: Lxx = (x*x).operator()
1234 sage: Qx = x.quadratic_representation()
1235 sage: Qy = y.quadratic_representation()
1236 sage: Qxy = x.quadratic_representation(y)
1237 sage: Qex = J.one().quadratic_representation(x)
1238 sage: n = ZZ.random_element(10)
1239 sage: Qxn = (x^n).quadratic_representation()
1243 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1246 Property 2 (multiply on the right for :trac:`28272`):
1248 sage: alpha = J.base_ring().random_element()
1249 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1254 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1257 sage: not x.is_invertible() or (
1260 ....: x.inverse().quadratic_representation() )
1263 sage: Qxy(J.one()) == x*y
1268 sage: not x.is_invertible() or (
1269 ....: x.quadratic_representation(x.inverse())*Qx
1270 ....: == Qx*x.quadratic_representation(x.inverse()) )
1273 sage: not x.is_invertible() or (
1274 ....: x.quadratic_representation(x.inverse())*Qx
1276 ....: 2*Lx*Qex - Qx )
1279 sage: 2*Lx*Qex - Qx == Lxx
1284 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1294 sage: not x.is_invertible() or (
1295 ....: Qx*x.inverse().operator() == Lx )
1300 sage: not x.operator_commutes_with(y) or (
1301 ....: Qx(y)^n == Qxn(y^n) )
1307 elif not other
in self
.parent():
1308 raise TypeError("'other' must live in the same algebra")
1311 M
= other
.operator()
1312 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1316 def spectral_decomposition(self
):
1318 Return the unique spectral decomposition of this element.
1322 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1323 element's left-multiplication-by operator to the subalgebra it
1324 generates. We then compute the spectral decomposition of that
1325 operator, and the spectral projectors we get back must be the
1326 left-multiplication-by operators for the idempotents we
1327 seek. Thus applying them to the identity element gives us those
1330 Since the eigenvalues are required to be distinct, we take
1331 the spectral decomposition of the zero element to be zero
1332 times the identity element of the algebra (which is idempotent,
1337 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1341 The spectral decomposition of the identity is ``1`` times itself,
1342 and the spectral decomposition of zero is ``0`` times the identity::
1344 sage: J = RealSymmetricEJA(3)
1347 sage: J.one().spectral_decomposition()
1349 sage: J.zero().spectral_decomposition()
1354 sage: J = RealSymmetricEJA(4)
1355 sage: x = sum(J.gens())
1356 sage: sd = x.spectral_decomposition()
1361 sage: c0.inner_product(c1) == 0
1363 sage: c0.is_idempotent()
1365 sage: c1.is_idempotent()
1367 sage: c0 + c1 == J.one()
1369 sage: l0*c0 + l1*c1 == x
1372 The spectral decomposition should work in subalgebras, too::
1374 sage: J = RealSymmetricEJA(4)
1375 sage: (e0, e1, e2, e3, e4, e5, e6, e7, e8, e9) = J.gens()
1376 sage: A = 2*e5 - 2*e8
1377 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1378 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1379 sage: (f0, f1, f2) = J1.gens()
1380 sage: f0.spectral_decomposition()
1384 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1386 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1387 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1390 def subalgebra_generated_by(self
, **kwargs
):
1392 Return the associative subalgebra of the parent EJA generated
1395 Since our parent algebra is unital, we want "subalgebra" to mean
1396 "unital subalgebra" as well; thus the subalgebra that an element
1397 generates will itself be a Euclidean Jordan algebra after
1398 restricting the algebra operations appropriately. This is the
1399 subalgebra that Faraut and Korányi work with in section II.2, for
1404 sage: from mjo.eja.eja_algebra import random_eja
1408 This subalgebra, being composed of only powers, is associative::
1410 sage: set_random_seed()
1411 sage: x0 = random_eja().random_element()
1412 sage: A = x0.subalgebra_generated_by()
1413 sage: x,y,z = A.random_elements(3)
1414 sage: (x*y)*z == x*(y*z)
1417 Squaring in the subalgebra should work the same as in
1420 sage: set_random_seed()
1421 sage: x = random_eja().random_element()
1422 sage: A = x.subalgebra_generated_by()
1423 sage: A(x^2) == A(x)*A(x)
1426 By definition, the subalgebra generated by the zero element is
1427 the one-dimensional algebra generated by the identity
1428 element... unless the original algebra was trivial, in which
1429 case the subalgebra is trivial too::
1431 sage: set_random_seed()
1432 sage: A = random_eja().zero().subalgebra_generated_by()
1433 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1437 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1438 A
= self
.parent().subalgebra(powers
, associative
=True, **kwargs
)
1439 A
.one
.set_cache(A(self
.parent().one()))
1443 def subalgebra_idempotent(self
):
1445 Find an idempotent in the associative subalgebra I generate
1446 using Proposition 2.3.5 in Baes.
1450 sage: from mjo.eja.eja_algebra import random_eja
1454 Ensure that we can find an idempotent in a non-trivial algebra
1455 where there are non-nilpotent elements, or that we get the dumb
1456 solution in the trivial algebra::
1458 sage: set_random_seed()
1459 sage: J = random_eja()
1460 sage: x = J.random_element()
1461 sage: while x.is_nilpotent() and not J.is_trivial():
1462 ....: x = J.random_element()
1463 sage: c = x.subalgebra_idempotent()
1468 if self
.parent().is_trivial():
1471 if self
.is_nilpotent():
1472 raise ValueError("this only works with non-nilpotent elements!")
1474 J
= self
.subalgebra_generated_by()
1477 # The image of the matrix of left-u^m-multiplication
1478 # will be minimal for some natural number s...
1480 minimal_dim
= J
.dimension()
1481 for i
in range(1, minimal_dim
):
1482 this_dim
= (u
**i
).operator().matrix().image().dimension()
1483 if this_dim
< minimal_dim
:
1484 minimal_dim
= this_dim
1487 # Now minimal_matrix should correspond to the smallest
1488 # non-zero subspace in Baes's (or really, Koecher's)
1491 # However, we need to restrict the matrix to work on the
1492 # subspace... or do we? Can't we just solve, knowing that
1493 # A(c) = u^(s+1) should have a solution in the big space,
1496 # Beware, solve_right() means that we're using COLUMN vectors.
1497 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1499 A
= u_next
.operator().matrix()
1500 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1502 # Now c is the idempotent we want, but it still lives in the subalgebra.
1503 return c
.superalgebra_element()
1508 Return my trace, the sum of my eigenvalues.
1510 In a trivial algebra, however you want to look at it, the trace is
1511 an empty sum for which we declare the result to be zero.
1515 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1522 sage: J = TrivialEJA()
1523 sage: J.zero().trace()
1527 sage: J = JordanSpinEJA(3)
1528 sage: x = sum(J.gens())
1534 sage: J = HadamardEJA(5)
1535 sage: J.one().trace()
1540 The trace of an element is a real number::
1542 sage: set_random_seed()
1543 sage: J = random_eja()
1544 sage: J.random_element().trace() in RLF
1547 The trace is linear::
1549 sage: set_random_seed()
1550 sage: J = random_eja()
1551 sage: x,y = J.random_elements(2)
1552 sage: alpha = J.base_ring().random_element()
1553 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1561 # Special case for the trivial algebra where
1562 # the trace is an empty sum.
1563 return P
.base_ring().zero()
1565 p
= P
._charpoly
_coefficients
()[r
-1]
1566 # The _charpoly_coeff function already adds the factor of
1567 # -1 to ensure that _charpoly_coeff(r-1) is really what
1568 # appears in front of t^{r-1} in the charpoly. However,
1569 # we want the negative of THAT for the trace.
1570 return -p(*self
.to_vector())
1573 def trace_inner_product(self
, other
):
1575 Return the trace inner product of myself and ``other``.
1579 sage: from mjo.eja.eja_algebra import random_eja
1583 The trace inner product is commutative, bilinear, and associative::
1585 sage: set_random_seed()
1586 sage: J = random_eja()
1587 sage: x,y,z = J.random_elements(3)
1589 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1592 sage: a = J.base_ring().random_element();
1593 sage: actual = (a*(x+z)).trace_inner_product(y)
1594 sage: expected = ( a*x.trace_inner_product(y) +
1595 ....: a*z.trace_inner_product(y) )
1596 sage: actual == expected
1598 sage: actual = x.trace_inner_product(a*(y+z))
1599 sage: expected = ( a*x.trace_inner_product(y) +
1600 ....: a*x.trace_inner_product(z) )
1601 sage: actual == expected
1604 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1608 if not other
in self
.parent():
1609 raise TypeError("'other' must live in the same algebra")
1611 return (self
*other
).trace()
1614 def trace_norm(self
):
1616 The norm of this element with respect to :meth:`trace_inner_product`.
1620 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1625 sage: J = HadamardEJA(2)
1626 sage: x = sum(J.gens())
1627 sage: x.trace_norm()
1632 sage: J = JordanSpinEJA(4)
1633 sage: x = sum(J.gens())
1634 sage: x.trace_norm()
1638 return self
.trace_inner_product(self
).sqrt()