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()) )
1123 The norm of this element with respect to :meth:`inner_product`.
1127 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1132 sage: J = HadamardEJA(2)
1133 sage: x = sum(J.gens())
1139 sage: J = JordanSpinEJA(4)
1140 sage: x = sum(J.gens())
1145 return self
.inner_product(self
).sqrt()
1150 Return the left-multiplication-by-this-element
1151 operator on the ambient algebra.
1155 sage: from mjo.eja.eja_algebra import random_eja
1159 sage: set_random_seed()
1160 sage: J = random_eja()
1161 sage: x,y = J.random_elements(2)
1162 sage: x.operator()(y) == x*y
1164 sage: y.operator()(x) == x*y
1169 left_mult_by_self
= lambda y
: self
*y
1170 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1171 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1174 def quadratic_representation(self
, other
=None):
1176 Return the quadratic representation of this element.
1180 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1185 The explicit form in the spin factor algebra is given by
1186 Alizadeh's Example 11.12::
1188 sage: set_random_seed()
1189 sage: x = JordanSpinEJA.random_instance().random_element()
1190 sage: x_vec = x.to_vector()
1191 sage: Q = matrix.identity(x.base_ring(), 0)
1192 sage: n = x_vec.degree()
1195 ....: x_bar = x_vec[1:]
1196 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1197 ....: B = 2*x0*x_bar.row()
1198 ....: C = 2*x0*x_bar.column()
1199 ....: D = matrix.identity(x.base_ring(), n-1)
1200 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1201 ....: D = D + 2*x_bar.tensor_product(x_bar)
1202 ....: Q = matrix.block(2,2,[A,B,C,D])
1203 sage: Q == x.quadratic_representation().matrix()
1206 Test all of the properties from Theorem 11.2 in Alizadeh::
1208 sage: set_random_seed()
1209 sage: J = random_eja()
1210 sage: x,y = J.random_elements(2)
1211 sage: Lx = x.operator()
1212 sage: Lxx = (x*x).operator()
1213 sage: Qx = x.quadratic_representation()
1214 sage: Qy = y.quadratic_representation()
1215 sage: Qxy = x.quadratic_representation(y)
1216 sage: Qex = J.one().quadratic_representation(x)
1217 sage: n = ZZ.random_element(10)
1218 sage: Qxn = (x^n).quadratic_representation()
1222 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1225 Property 2 (multiply on the right for :trac:`28272`):
1227 sage: alpha = J.base_ring().random_element()
1228 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1233 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1236 sage: not x.is_invertible() or (
1239 ....: x.inverse().quadratic_representation() )
1242 sage: Qxy(J.one()) == x*y
1247 sage: not x.is_invertible() or (
1248 ....: x.quadratic_representation(x.inverse())*Qx
1249 ....: == Qx*x.quadratic_representation(x.inverse()) )
1252 sage: not x.is_invertible() or (
1253 ....: x.quadratic_representation(x.inverse())*Qx
1255 ....: 2*Lx*Qex - Qx )
1258 sage: 2*Lx*Qex - Qx == Lxx
1263 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1273 sage: not x.is_invertible() or (
1274 ....: Qx*x.inverse().operator() == Lx )
1279 sage: not x.operator_commutes_with(y) or (
1280 ....: Qx(y)^n == Qxn(y^n) )
1286 elif not other
in self
.parent():
1287 raise TypeError("'other' must live in the same algebra")
1290 M
= other
.operator()
1291 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1295 def spectral_decomposition(self
):
1297 Return the unique spectral decomposition of this element.
1301 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1302 element's left-multiplication-by operator to the subalgebra it
1303 generates. We then compute the spectral decomposition of that
1304 operator, and the spectral projectors we get back must be the
1305 left-multiplication-by operators for the idempotents we
1306 seek. Thus applying them to the identity element gives us those
1309 Since the eigenvalues are required to be distinct, we take
1310 the spectral decomposition of the zero element to be zero
1311 times the identity element of the algebra (which is idempotent,
1316 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1320 The spectral decomposition of the identity is ``1`` times itself,
1321 and the spectral decomposition of zero is ``0`` times the identity::
1323 sage: J = RealSymmetricEJA(3)
1326 sage: J.one().spectral_decomposition()
1328 sage: J.zero().spectral_decomposition()
1333 sage: J = RealSymmetricEJA(4)
1334 sage: x = sum(J.gens())
1335 sage: sd = x.spectral_decomposition()
1340 sage: c0.inner_product(c1) == 0
1342 sage: c0.is_idempotent()
1344 sage: c1.is_idempotent()
1346 sage: c0 + c1 == J.one()
1348 sage: l0*c0 + l1*c1 == x
1351 The spectral decomposition should work in subalgebras, too::
1353 sage: J = RealSymmetricEJA(4)
1354 sage: (e0, e1, e2, e3, e4, e5, e6, e7, e8, e9) = J.gens()
1355 sage: A = 2*e5 - 2*e8
1356 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1357 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1358 sage: (f0, f1, f2) = J1.gens()
1359 sage: f0.spectral_decomposition()
1363 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1365 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1366 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1369 def subalgebra_generated_by(self
, **kwargs
):
1371 Return the associative subalgebra of the parent EJA generated
1374 Since our parent algebra is unital, we want "subalgebra" to mean
1375 "unital subalgebra" as well; thus the subalgebra that an element
1376 generates will itself be a Euclidean Jordan algebra after
1377 restricting the algebra operations appropriately. This is the
1378 subalgebra that Faraut and Korányi work with in section II.2, for
1383 sage: from mjo.eja.eja_algebra import random_eja
1387 This subalgebra, being composed of only powers, is associative::
1389 sage: set_random_seed()
1390 sage: x0 = random_eja().random_element()
1391 sage: A = x0.subalgebra_generated_by()
1392 sage: x,y,z = A.random_elements(3)
1393 sage: (x*y)*z == x*(y*z)
1396 Squaring in the subalgebra should work the same as in
1399 sage: set_random_seed()
1400 sage: x = random_eja().random_element()
1401 sage: A = x.subalgebra_generated_by()
1402 sage: A(x^2) == A(x)*A(x)
1405 By definition, the subalgebra generated by the zero element is
1406 the one-dimensional algebra generated by the identity
1407 element... unless the original algebra was trivial, in which
1408 case the subalgebra is trivial too::
1410 sage: set_random_seed()
1411 sage: A = random_eja().zero().subalgebra_generated_by()
1412 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1416 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1417 A
= self
.parent().subalgebra(powers
, associative
=True, **kwargs
)
1418 A
.one
.set_cache(A(self
.parent().one()))
1422 def subalgebra_idempotent(self
):
1424 Find an idempotent in the associative subalgebra I generate
1425 using Proposition 2.3.5 in Baes.
1429 sage: from mjo.eja.eja_algebra import random_eja
1433 Ensure that we can find an idempotent in a non-trivial algebra
1434 where there are non-nilpotent elements, or that we get the dumb
1435 solution in the trivial algebra::
1437 sage: set_random_seed()
1438 sage: J = random_eja()
1439 sage: x = J.random_element()
1440 sage: while x.is_nilpotent() and not J.is_trivial():
1441 ....: x = J.random_element()
1442 sage: c = x.subalgebra_idempotent()
1447 if self
.parent().is_trivial():
1450 if self
.is_nilpotent():
1451 raise ValueError("this only works with non-nilpotent elements!")
1453 J
= self
.subalgebra_generated_by()
1456 # The image of the matrix of left-u^m-multiplication
1457 # will be minimal for some natural number s...
1459 minimal_dim
= J
.dimension()
1460 for i
in range(1, minimal_dim
):
1461 this_dim
= (u
**i
).operator().matrix().image().dimension()
1462 if this_dim
< minimal_dim
:
1463 minimal_dim
= this_dim
1466 # Now minimal_matrix should correspond to the smallest
1467 # non-zero subspace in Baes's (or really, Koecher's)
1470 # However, we need to restrict the matrix to work on the
1471 # subspace... or do we? Can't we just solve, knowing that
1472 # A(c) = u^(s+1) should have a solution in the big space,
1475 # Beware, solve_right() means that we're using COLUMN vectors.
1476 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1478 A
= u_next
.operator().matrix()
1479 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1481 # Now c is the idempotent we want, but it still lives in the subalgebra.
1482 return c
.superalgebra_element()
1487 Return my trace, the sum of my eigenvalues.
1489 In a trivial algebra, however you want to look at it, the trace is
1490 an empty sum for which we declare the result to be zero.
1494 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1501 sage: J = TrivialEJA()
1502 sage: J.zero().trace()
1506 sage: J = JordanSpinEJA(3)
1507 sage: x = sum(J.gens())
1513 sage: J = HadamardEJA(5)
1514 sage: J.one().trace()
1519 The trace of an element is a real number::
1521 sage: set_random_seed()
1522 sage: J = random_eja()
1523 sage: J.random_element().trace() in RLF
1526 The trace is linear::
1528 sage: set_random_seed()
1529 sage: J = random_eja()
1530 sage: x,y = J.random_elements(2)
1531 sage: alpha = J.base_ring().random_element()
1532 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1540 # Special case for the trivial algebra where
1541 # the trace is an empty sum.
1542 return P
.base_ring().zero()
1544 p
= P
._charpoly
_coefficients
()[r
-1]
1545 # The _charpoly_coeff function already adds the factor of
1546 # -1 to ensure that _charpoly_coeff(r-1) is really what
1547 # appears in front of t^{r-1} in the charpoly. However,
1548 # we want the negative of THAT for the trace.
1549 return -p(*self
.to_vector())
1552 def trace_inner_product(self
, other
):
1554 Return the trace inner product of myself and ``other``.
1558 sage: from mjo.eja.eja_algebra import random_eja
1562 The trace inner product is commutative, bilinear, and associative::
1564 sage: set_random_seed()
1565 sage: J = random_eja()
1566 sage: x,y,z = J.random_elements(3)
1568 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1571 sage: a = J.base_ring().random_element();
1572 sage: actual = (a*(x+z)).trace_inner_product(y)
1573 sage: expected = ( a*x.trace_inner_product(y) +
1574 ....: a*z.trace_inner_product(y) )
1575 sage: actual == expected
1577 sage: actual = x.trace_inner_product(a*(y+z))
1578 sage: expected = ( a*x.trace_inner_product(y) +
1579 ....: a*x.trace_inner_product(z) )
1580 sage: actual == expected
1583 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1587 if not other
in self
.parent():
1588 raise TypeError("'other' must live in the same algebra")
1590 return (self
*other
).trace()
1593 def trace_norm(self
):
1595 The norm of this element with respect to :meth:`trace_inner_product`.
1599 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1604 sage: J = HadamardEJA(2)
1605 sage: x = sum(J.gens())
1606 sage: x.trace_norm()
1611 sage: J = JordanSpinEJA(4)
1612 sage: x = sum(J.gens())
1613 sage: x.trace_norm()
1617 return self
.trace_inner_product(self
).sqrt()
1621 class CartesianProductEJAElement(FiniteDimensionalEJAElement
):
1623 def to_matrix(self
):
1627 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1628 ....: RealSymmetricEJA)
1632 sage: J1 = HadamardEJA(1)
1633 sage: J2 = RealSymmetricEJA(2)
1634 sage: J = cartesian_product([J1,J2])
1635 sage: x = sum(J.gens())
1636 sage: x.to_matrix()[0]
1638 sage: x.to_matrix()[1]
1639 [ 1 0.7071067811865475?]
1640 [0.7071067811865475? 1]
1643 B
= self
.parent().matrix_basis()
1644 W
= self
.parent().matrix_space()
1646 # Aaaaand linear combinations don't work in Cartesian
1647 # product spaces, even though they provide a method
1649 pairs
= zip(B
, self
.to_vector())
1650 return sum( ( W(tuple(alpha
*b_i
for b_i
in b
))
1651 for (b
,alpha
) in pairs
),