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
, _scale
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.monomial(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,
1081 ....: QuaternionHermitianEJA,
1082 ....: RealSymmetricEJA)
1086 sage: J = ComplexHermitianEJA(3)
1089 sage: J.one().to_matrix()
1099 sage: J = QuaternionHermitianEJA(2)
1102 sage: J.one().to_matrix()
1112 This also works in Cartesian product algebras::
1114 sage: J1 = HadamardEJA(1)
1115 sage: J2 = RealSymmetricEJA(2)
1116 sage: J = cartesian_product([J1,J2])
1117 sage: x = sum(J.gens())
1118 sage: x.to_matrix()[0]
1120 sage: x.to_matrix()[1]
1121 [ 1 0.7071067811865475?]
1122 [0.7071067811865475? 1]
1125 B
= self
.parent().matrix_basis()
1126 W
= self
.parent().matrix_space()
1128 if self
.parent()._matrix
_basis
_is
_cartesian
:
1129 # Aaaaand linear combinations don't work in Cartesian
1130 # product spaces, even though they provide a method
1131 # with that name. This is special-cased because the
1132 # _scale() function is slow.
1133 pairs
= zip(B
, self
.to_vector())
1134 return sum( ( _scale(b
, alpha
) for (b
,alpha
) in pairs
),
1137 # This is just a manual "from_vector()", but of course
1138 # matrix spaces aren't vector spaces in sage, so they
1139 # don't have a from_vector() method.
1140 return W
.linear_combination( zip(B
, self
.to_vector()) )
1146 The norm of this element with respect to :meth:`inner_product`.
1150 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1155 sage: J = HadamardEJA(2)
1156 sage: x = sum(J.gens())
1162 sage: J = JordanSpinEJA(4)
1163 sage: x = sum(J.gens())
1168 return self
.inner_product(self
).sqrt()
1173 Return the left-multiplication-by-this-element
1174 operator on the ambient algebra.
1178 sage: from mjo.eja.eja_algebra import random_eja
1182 sage: set_random_seed()
1183 sage: J = random_eja()
1184 sage: x,y = J.random_elements(2)
1185 sage: x.operator()(y) == x*y
1187 sage: y.operator()(x) == x*y
1192 left_mult_by_self
= lambda y
: self
*y
1193 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1194 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1197 def quadratic_representation(self
, other
=None):
1199 Return the quadratic representation of this element.
1203 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1208 The explicit form in the spin factor algebra is given by
1209 Alizadeh's Example 11.12::
1211 sage: set_random_seed()
1212 sage: x = JordanSpinEJA.random_instance().random_element()
1213 sage: x_vec = x.to_vector()
1214 sage: Q = matrix.identity(x.base_ring(), 0)
1215 sage: n = x_vec.degree()
1218 ....: x_bar = x_vec[1:]
1219 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1220 ....: B = 2*x0*x_bar.row()
1221 ....: C = 2*x0*x_bar.column()
1222 ....: D = matrix.identity(x.base_ring(), n-1)
1223 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1224 ....: D = D + 2*x_bar.tensor_product(x_bar)
1225 ....: Q = matrix.block(2,2,[A,B,C,D])
1226 sage: Q == x.quadratic_representation().matrix()
1229 Test all of the properties from Theorem 11.2 in Alizadeh::
1231 sage: set_random_seed()
1232 sage: J = random_eja()
1233 sage: x,y = J.random_elements(2)
1234 sage: Lx = x.operator()
1235 sage: Lxx = (x*x).operator()
1236 sage: Qx = x.quadratic_representation()
1237 sage: Qy = y.quadratic_representation()
1238 sage: Qxy = x.quadratic_representation(y)
1239 sage: Qex = J.one().quadratic_representation(x)
1240 sage: n = ZZ.random_element(10)
1241 sage: Qxn = (x^n).quadratic_representation()
1245 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1248 Property 2 (multiply on the right for :trac:`28272`):
1250 sage: alpha = J.base_ring().random_element()
1251 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1256 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1259 sage: not x.is_invertible() or (
1262 ....: x.inverse().quadratic_representation() )
1265 sage: Qxy(J.one()) == x*y
1270 sage: not x.is_invertible() or (
1271 ....: x.quadratic_representation(x.inverse())*Qx
1272 ....: == Qx*x.quadratic_representation(x.inverse()) )
1275 sage: not x.is_invertible() or (
1276 ....: x.quadratic_representation(x.inverse())*Qx
1278 ....: 2*Lx*Qex - Qx )
1281 sage: 2*Lx*Qex - Qx == Lxx
1286 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1296 sage: not x.is_invertible() or (
1297 ....: Qx*x.inverse().operator() == Lx )
1302 sage: not x.operator_commutes_with(y) or (
1303 ....: Qx(y)^n == Qxn(y^n) )
1309 elif not other
in self
.parent():
1310 raise TypeError("'other' must live in the same algebra")
1313 M
= other
.operator()
1314 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1318 def spectral_decomposition(self
):
1320 Return the unique spectral decomposition of this element.
1324 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1325 element's left-multiplication-by operator to the subalgebra it
1326 generates. We then compute the spectral decomposition of that
1327 operator, and the spectral projectors we get back must be the
1328 left-multiplication-by operators for the idempotents we
1329 seek. Thus applying them to the identity element gives us those
1332 Since the eigenvalues are required to be distinct, we take
1333 the spectral decomposition of the zero element to be zero
1334 times the identity element of the algebra (which is idempotent,
1339 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1343 The spectral decomposition of the identity is ``1`` times itself,
1344 and the spectral decomposition of zero is ``0`` times the identity::
1346 sage: J = RealSymmetricEJA(3)
1349 sage: J.one().spectral_decomposition()
1351 sage: J.zero().spectral_decomposition()
1356 sage: J = RealSymmetricEJA(4)
1357 sage: x = sum(J.gens())
1358 sage: sd = x.spectral_decomposition()
1363 sage: c0.inner_product(c1) == 0
1365 sage: c0.is_idempotent()
1367 sage: c1.is_idempotent()
1369 sage: c0 + c1 == J.one()
1371 sage: l0*c0 + l1*c1 == x
1374 The spectral decomposition should work in subalgebras, too::
1376 sage: J = RealSymmetricEJA(4)
1377 sage: (e0, e1, e2, e3, e4, e5, e6, e7, e8, e9) = J.gens()
1378 sage: A = 2*e5 - 2*e8
1379 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1380 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1381 sage: (f0, f1, f2) = J1.gens()
1382 sage: f0.spectral_decomposition()
1386 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1388 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1389 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1392 def subalgebra_generated_by(self
, **kwargs
):
1394 Return the associative subalgebra of the parent EJA generated
1397 Since our parent algebra is unital, we want "subalgebra" to mean
1398 "unital subalgebra" as well; thus the subalgebra that an element
1399 generates will itself be a Euclidean Jordan algebra after
1400 restricting the algebra operations appropriately. This is the
1401 subalgebra that Faraut and Korányi work with in section II.2, for
1406 sage: from mjo.eja.eja_algebra import (random_eja,
1408 ....: RealSymmetricEJA)
1412 We can create subalgebras of Cartesian product EJAs that are not
1413 themselves Cartesian product EJAs (they're just "regular" EJAs)::
1415 sage: J1 = HadamardEJA(3)
1416 sage: J2 = RealSymmetricEJA(2)
1417 sage: J = cartesian_product([J1,J2])
1418 sage: J.one().subalgebra_generated_by()
1419 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
1423 This subalgebra, being composed of only powers, is associative::
1425 sage: set_random_seed()
1426 sage: x0 = random_eja().random_element()
1427 sage: A = x0.subalgebra_generated_by()
1428 sage: x,y,z = A.random_elements(3)
1429 sage: (x*y)*z == x*(y*z)
1432 Squaring in the subalgebra should work the same as in
1435 sage: set_random_seed()
1436 sage: x = random_eja().random_element()
1437 sage: A = x.subalgebra_generated_by()
1438 sage: A(x^2) == A(x)*A(x)
1441 By definition, the subalgebra generated by the zero element is
1442 the one-dimensional algebra generated by the identity
1443 element... unless the original algebra was trivial, in which
1444 case the subalgebra is trivial too::
1446 sage: set_random_seed()
1447 sage: A = random_eja().zero().subalgebra_generated_by()
1448 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1452 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1453 A
= self
.parent().subalgebra(powers
,
1458 A
.one
.set_cache(A(self
.parent().one()))
1462 def subalgebra_idempotent(self
):
1464 Find an idempotent in the associative subalgebra I generate
1465 using Proposition 2.3.5 in Baes.
1469 sage: from mjo.eja.eja_algebra import random_eja
1473 Ensure that we can find an idempotent in a non-trivial algebra
1474 where there are non-nilpotent elements, or that we get the dumb
1475 solution in the trivial algebra::
1477 sage: set_random_seed()
1478 sage: J = random_eja()
1479 sage: x = J.random_element()
1480 sage: while x.is_nilpotent() and not J.is_trivial():
1481 ....: x = J.random_element()
1482 sage: c = x.subalgebra_idempotent()
1487 if self
.parent().is_trivial():
1490 if self
.is_nilpotent():
1491 raise ValueError("this only works with non-nilpotent elements!")
1493 J
= self
.subalgebra_generated_by()
1496 # The image of the matrix of left-u^m-multiplication
1497 # will be minimal for some natural number s...
1499 minimal_dim
= J
.dimension()
1500 for i
in range(1, minimal_dim
):
1501 this_dim
= (u
**i
).operator().matrix().image().dimension()
1502 if this_dim
< minimal_dim
:
1503 minimal_dim
= this_dim
1506 # Now minimal_matrix should correspond to the smallest
1507 # non-zero subspace in Baes's (or really, Koecher's)
1510 # However, we need to restrict the matrix to work on the
1511 # subspace... or do we? Can't we just solve, knowing that
1512 # A(c) = u^(s+1) should have a solution in the big space,
1515 # Beware, solve_right() means that we're using COLUMN vectors.
1516 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1518 A
= u_next
.operator().matrix()
1519 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1521 # Now c is the idempotent we want, but it still lives in the subalgebra.
1522 return c
.superalgebra_element()
1527 Return my trace, the sum of my eigenvalues.
1529 In a trivial algebra, however you want to look at it, the trace is
1530 an empty sum for which we declare the result to be zero.
1534 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1541 sage: J = TrivialEJA()
1542 sage: J.zero().trace()
1546 sage: J = JordanSpinEJA(3)
1547 sage: x = sum(J.gens())
1553 sage: J = HadamardEJA(5)
1554 sage: J.one().trace()
1559 The trace of an element is a real number::
1561 sage: set_random_seed()
1562 sage: J = random_eja()
1563 sage: J.random_element().trace() in RLF
1566 The trace is linear::
1568 sage: set_random_seed()
1569 sage: J = random_eja()
1570 sage: x,y = J.random_elements(2)
1571 sage: alpha = J.base_ring().random_element()
1572 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1580 # Special case for the trivial algebra where
1581 # the trace is an empty sum.
1582 return P
.base_ring().zero()
1584 p
= P
._charpoly
_coefficients
()[r
-1]
1585 # The _charpoly_coeff function already adds the factor of
1586 # -1 to ensure that _charpoly_coeff(r-1) is really what
1587 # appears in front of t^{r-1} in the charpoly. However,
1588 # we want the negative of THAT for the trace.
1589 return -p(*self
.to_vector())
1592 def trace_inner_product(self
, other
):
1594 Return the trace inner product of myself and ``other``.
1598 sage: from mjo.eja.eja_algebra import random_eja
1602 The trace inner product is commutative, bilinear, and associative::
1604 sage: set_random_seed()
1605 sage: J = random_eja()
1606 sage: x,y,z = J.random_elements(3)
1608 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1611 sage: a = J.base_ring().random_element();
1612 sage: actual = (a*(x+z)).trace_inner_product(y)
1613 sage: expected = ( a*x.trace_inner_product(y) +
1614 ....: a*z.trace_inner_product(y) )
1615 sage: actual == expected
1617 sage: actual = x.trace_inner_product(a*(y+z))
1618 sage: expected = ( a*x.trace_inner_product(y) +
1619 ....: a*x.trace_inner_product(z) )
1620 sage: actual == expected
1623 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1627 if not other
in self
.parent():
1628 raise TypeError("'other' must live in the same algebra")
1630 return (self
*other
).trace()
1633 def trace_norm(self
):
1635 The norm of this element with respect to :meth:`trace_inner_product`.
1639 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1644 sage: J = HadamardEJA(2)
1645 sage: x = sum(J.gens())
1646 sage: x.trace_norm()
1651 sage: J = JordanSpinEJA(4)
1652 sage: x = sum(J.gens())
1653 sage: x.trace_norm()
1657 return self
.trace_inner_product(self
).sqrt()