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: b0,b1,b2 = J.gens()
171 sage: A = (b0 + 2*b1 + 3*b2).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: x = sum( J.gens() )
357 sage: J = JordanSpinEJA(3)
358 sage: x = sum( J.gens() )
362 The determinant of the sole element in the rank-zero trivial
363 algebra is ``1``, by three paths of reasoning. First, its
364 characteristic polynomial is a constant ``1``, so the constant
365 term in that polynomial is ``1``. Second, the characteristic
366 polynomial evaluated at zero is again ``1``. And finally, the
367 (empty) product of its eigenvalues is likewise just unity::
369 sage: J = TrivialEJA()
375 An element is invertible if and only if its determinant is
378 sage: set_random_seed()
379 sage: x = random_eja().random_element()
380 sage: x.is_invertible() == (x.det() != 0)
383 Ensure that the determinant is multiplicative on an associative
384 subalgebra as in Faraut and Korányi's Proposition II.2.2::
386 sage: set_random_seed()
387 sage: J = random_eja().random_element().subalgebra_generated_by()
388 sage: x,y = J.random_elements(2)
389 sage: (x*y).det() == x.det()*y.det()
392 The determinant in matrix algebras is just the usual determinant::
394 sage: set_random_seed()
395 sage: X = matrix.random(QQ,3)
397 sage: J1 = RealSymmetricEJA(3)
398 sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
399 sage: expected = X.det()
400 sage: actual1 = J1(X).det()
401 sage: actual2 = J2(X).det()
402 sage: actual1 == expected
404 sage: actual2 == expected
409 sage: set_random_seed()
410 sage: J1 = ComplexHermitianEJA(2)
411 sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
412 sage: X = matrix.random(GaussianIntegers(), 2)
414 sage: expected = AA(X.det())
415 sage: actual1 = J1(J1.real_embed(X)).det()
416 sage: actual2 = J2(J2.real_embed(X)).det()
417 sage: expected == actual1
419 sage: expected == actual2
427 # Special case, since we don't get the a0=1
428 # coefficient when the rank of the algebra
430 return P
.base_ring().one()
432 p
= P
._charpoly
_coefficients
()[0]
433 # The _charpoly_coeff function already adds the factor of -1
434 # to ensure that _charpoly_coefficients()[0] is really what
435 # appears in front of t^{0} in the charpoly. However, we want
436 # (-1)^r times THAT for the determinant.
437 return ((-1)**r
)*p(*self
.to_vector())
443 Return the Jordan-multiplicative inverse of this element.
447 In general we appeal to the quadratic representation as in
448 Koecher's Theorem 12 in Chapter III, Section 5. But if the
449 parent algebra's "characteristic polynomial of" coefficients
450 happen to be cached, then we use Proposition II.2.4 in Faraut
451 and Korányi which gives a formula for the inverse based on the
452 characteristic polynomial and the Cayley-Hamilton theorem for
453 Euclidean Jordan algebras::
457 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
463 The inverse in the spin factor algebra is given in Alizadeh's
466 sage: set_random_seed()
467 sage: J = JordanSpinEJA.random_instance()
468 sage: x = J.random_element()
469 sage: while not x.is_invertible():
470 ....: x = J.random_element()
471 sage: x_vec = x.to_vector()
473 sage: x_bar = x_vec[1:]
474 sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
475 sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
476 sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
477 sage: x.inverse() == J.from_vector(x_inverse)
480 Trying to invert a non-invertible element throws an error:
482 sage: JordanSpinEJA(3).zero().inverse()
483 Traceback (most recent call last):
485 ZeroDivisionError: element is not invertible
489 The identity element is its own inverse::
491 sage: set_random_seed()
492 sage: J = random_eja()
493 sage: J.one().inverse() == J.one()
496 If an element has an inverse, it acts like one::
498 sage: set_random_seed()
499 sage: J = random_eja()
500 sage: x = J.random_element()
501 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
504 The inverse of the inverse is what we started with::
506 sage: set_random_seed()
507 sage: J = random_eja()
508 sage: x = J.random_element()
509 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
512 Proposition II.2.3 in Faraut and Korányi says that the inverse
513 of an element is the inverse of its left-multiplication operator
514 applied to the algebra's identity, when that inverse exists::
516 sage: set_random_seed()
517 sage: J = random_eja()
518 sage: x = J.random_element()
519 sage: (not x.operator().is_invertible()) or (
520 ....: x.operator().inverse()(J.one()) == x.inverse() )
523 Check that the fast (cached) and slow algorithms give the same
526 sage: set_random_seed() # long time
527 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
528 sage: x = J.random_element() # long time
529 sage: while not x.is_invertible(): # long time
530 ....: x = J.random_element() # long time
531 sage: slow = x.inverse() # long time
532 sage: _ = J._charpoly_coefficients() # long time
533 sage: fast = x.inverse() # long time
534 sage: slow == fast # long time
537 not_invertible_msg
= "element is not invertible"
538 if self
.parent()._charpoly
_coefficients
.is_in_cache():
539 # We can invert using our charpoly if it will be fast to
540 # compute. If the coefficients are cached, our rank had
542 if self
.det().is_zero():
543 raise ZeroDivisionError(not_invertible_msg
)
544 r
= self
.parent().rank()
545 a
= self
.characteristic_polynomial().coefficients(sparse
=False)
546 return (-1)**(r
+1)*sum(a
[i
+1]*self
**i
for i
in range(r
))/self
.det()
549 inv
= (~self
.quadratic_representation())(self
)
550 self
.is_invertible
.set_cache(True)
552 except ZeroDivisionError:
553 self
.is_invertible
.set_cache(False)
554 raise ZeroDivisionError(not_invertible_msg
)
558 def is_invertible(self
):
560 Return whether or not this element is invertible.
564 If computing my determinant will be fast, we do so and compare
565 with zero (Proposition II.2.4 in Faraut and
566 Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
567 reduces the problem to the invertibility of my quadratic
572 sage: from mjo.eja.eja_algebra import random_eja
576 The identity element is always invertible::
578 sage: set_random_seed()
579 sage: J = random_eja()
580 sage: J.one().is_invertible()
583 The zero element is never invertible in a non-trivial algebra::
585 sage: set_random_seed()
586 sage: J = random_eja()
587 sage: (not J.is_trivial()) and J.zero().is_invertible()
590 Test that the fast (cached) and slow algorithms give the same
593 sage: set_random_seed() # long time
594 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
595 sage: x = J.random_element() # long time
596 sage: slow = x.is_invertible() # long time
597 sage: _ = J._charpoly_coefficients() # long time
598 sage: fast = x.is_invertible() # long time
599 sage: slow == fast # long time
603 if self
.parent().is_trivial():
608 if self
.parent()._charpoly
_coefficients
.is_in_cache():
609 # The determinant will be quicker than inverting the
610 # quadratic representation, most likely.
611 return (not self
.det().is_zero())
613 # The easiest way to determine if I'm invertible is to try.
615 inv
= (~self
.quadratic_representation())(self
)
616 self
.inverse
.set_cache(inv
)
618 except ZeroDivisionError:
622 def is_primitive_idempotent(self
):
624 Return whether or not this element is a primitive (or minimal)
627 A primitive idempotent is a non-zero idempotent that is not
628 the sum of two other non-zero idempotents. Remark 2.7.15 in
629 Baes shows that this is what he refers to as a "minimal
632 An element of a Euclidean Jordan algebra is a minimal idempotent
633 if it :meth:`is_idempotent` and if its Peirce subalgebra
634 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
639 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
640 ....: RealSymmetricEJA,
646 This method is sloooooow.
650 The spectral decomposition of a non-regular element should always
651 contain at least one non-minimal idempotent::
653 sage: J = RealSymmetricEJA(3)
654 sage: x = sum(J.gens())
657 sage: [ c.is_primitive_idempotent()
658 ....: for (l,c) in x.spectral_decomposition() ]
661 On the other hand, the spectral decomposition of a regular
662 element should always be in terms of minimal idempotents::
664 sage: J = JordanSpinEJA(4)
665 sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
668 sage: [ c.is_primitive_idempotent()
669 ....: for (l,c) in x.spectral_decomposition() ]
674 The identity element is minimal only in an EJA of rank one::
676 sage: set_random_seed()
677 sage: J = random_eja()
678 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
681 A non-idempotent cannot be a minimal idempotent::
683 sage: set_random_seed()
684 sage: J = JordanSpinEJA(4)
685 sage: x = J.random_element()
686 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
689 Proposition 2.7.19 in Baes says that an element is a minimal
690 idempotent if and only if it's idempotent with trace equal to
693 sage: set_random_seed()
694 sage: J = JordanSpinEJA(4)
695 sage: x = J.random_element()
696 sage: expected = (x.is_idempotent() and x.trace() == 1)
697 sage: actual = x.is_primitive_idempotent()
698 sage: actual == expected
701 Primitive idempotents must be non-zero::
703 sage: set_random_seed()
704 sage: J = random_eja()
705 sage: J.zero().is_idempotent()
707 sage: J.zero().is_primitive_idempotent()
710 As a consequence of the fact that primitive idempotents must
711 be non-zero, there are no primitive idempotents in a trivial
712 Euclidean Jordan algebra::
714 sage: J = TrivialEJA()
715 sage: J.one().is_idempotent()
717 sage: J.one().is_primitive_idempotent()
721 if not self
.is_idempotent():
727 (_
,_
,J1
) = self
.parent().peirce_decomposition(self
)
728 return (J1
.dimension() == 1)
731 def is_nilpotent(self
):
733 Return whether or not some power of this element is zero.
737 We use Theorem 5 in Chapter III of Koecher, which says that
738 an element ``x`` is nilpotent if and only if ``x.operator()``
739 is nilpotent. And it is a basic fact of linear algebra that
740 an operator on an `n`-dimensional space is nilpotent if and
741 only if, when raised to the `n`th power, it equals the zero
742 operator (for example, see Axler Corollary 8.8).
746 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
751 sage: J = JordanSpinEJA(3)
752 sage: x = sum(J.gens())
753 sage: x.is_nilpotent()
758 The identity element is never nilpotent, except in a trivial EJA::
760 sage: set_random_seed()
761 sage: J = random_eja()
762 sage: J.one().is_nilpotent() and not J.is_trivial()
765 The additive identity is always nilpotent::
767 sage: set_random_seed()
768 sage: random_eja().zero().is_nilpotent()
773 zero_operator
= P
.zero().operator()
774 return self
.operator()**P
.dimension() == zero_operator
777 def is_regular(self
):
779 Return whether or not this is a regular element.
783 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
788 The identity element always has degree one, but any element
789 linearly-independent from it is regular::
791 sage: J = JordanSpinEJA(5)
792 sage: J.one().is_regular()
794 sage: b0, b1, b2, b3, b4 = J.gens()
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: b0,b1,b2,b3 = J.gens()
847 sage: (b0 - b1).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 # Pretty sure we know what the minimal polynomial of
1051 # the zero operator is going to be. This ensures
1052 # consistency of e.g. the polynomial variable returned
1053 # in the "normal" case without us having to think about it.
1054 return self
.operator().minimal_polynomial()
1056 # If we don't orthonormalize the subalgebra's basis, then the
1057 # first two monomials in the subalgebra will be self^0 and
1058 # self^1... assuming that self^1 is not a scalar multiple of
1059 # self^0 (the unit element). We special case these to avoid
1060 # having to solve a system to coerce self into the subalgebra.
1061 A
= self
.subalgebra_generated_by(orthonormalize
=False)
1063 if A
.dimension() == 1:
1064 # Does a solve to find the scalar multiple alpha such that
1065 # alpha*unit = self. We have to do this because the basis
1066 # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
1067 unit
= self
.parent().one()
1068 alpha
= self
.to_vector() / unit
.to_vector()
1069 return (unit
.operator()*alpha
).minimal_polynomial()
1071 # If the dimension of the subalgebra is >= 2, then we just
1072 # use the second basis element.
1073 return A
.monomial(1).operator().minimal_polynomial()
1077 def to_matrix(self
):
1079 Return an (often more natural) representation of this element as a
1082 Every finite-dimensional Euclidean Jordan Algebra is a direct
1083 sum of five simple algebras, four of which comprise Hermitian
1084 matrices. This method returns a "natural" matrix
1085 representation of this element as either a Hermitian matrix or
1090 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1092 ....: QuaternionHermitianEJA,
1093 ....: RealSymmetricEJA)
1097 sage: J = ComplexHermitianEJA(3)
1100 sage: J.one().to_matrix()
1110 sage: J = QuaternionHermitianEJA(2)
1113 sage: J.one().to_matrix()
1123 This also works in Cartesian product algebras::
1125 sage: J1 = HadamardEJA(1)
1126 sage: J2 = RealSymmetricEJA(2)
1127 sage: J = cartesian_product([J1,J2])
1128 sage: x = sum(J.gens())
1129 sage: x.to_matrix()[0]
1131 sage: x.to_matrix()[1]
1132 [ 1 0.7071067811865475?]
1133 [0.7071067811865475? 1]
1136 B
= self
.parent().matrix_basis()
1137 W
= self
.parent().matrix_space()
1139 if hasattr(W
, 'cartesian_factors'):
1140 # Aaaaand linear combinations don't work in Cartesian
1141 # product spaces, even though they provide a method with
1142 # that name. This is hidden behind an "if" because the
1143 # _scale() function is slow.
1144 pairs
= zip(B
, self
.to_vector())
1145 return W
.sum( _scale(b
, alpha
) for (b
,alpha
) in pairs
)
1147 # This is just a manual "from_vector()", but of course
1148 # matrix spaces aren't vector spaces in sage, so they
1149 # don't have a from_vector() method.
1150 return W
.linear_combination( zip(B
, self
.to_vector()) )
1156 The norm of this element with respect to :meth:`inner_product`.
1160 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1165 sage: J = HadamardEJA(2)
1166 sage: x = sum(J.gens())
1172 sage: J = JordanSpinEJA(4)
1173 sage: x = sum(J.gens())
1178 return self
.inner_product(self
).sqrt()
1183 Return the left-multiplication-by-this-element
1184 operator on the ambient algebra.
1188 sage: from mjo.eja.eja_algebra import random_eja
1192 sage: set_random_seed()
1193 sage: J = random_eja()
1194 sage: x,y = J.random_elements(2)
1195 sage: x.operator()(y) == x*y
1197 sage: y.operator()(x) == x*y
1202 left_mult_by_self
= lambda y
: self
*y
1203 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1204 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1207 def quadratic_representation(self
, other
=None):
1209 Return the quadratic representation of this element.
1213 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1218 The explicit form in the spin factor algebra is given by
1219 Alizadeh's Example 11.12::
1221 sage: set_random_seed()
1222 sage: x = JordanSpinEJA.random_instance().random_element()
1223 sage: x_vec = x.to_vector()
1224 sage: Q = matrix.identity(x.base_ring(), 0)
1225 sage: n = x_vec.degree()
1228 ....: x_bar = x_vec[1:]
1229 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1230 ....: B = 2*x0*x_bar.row()
1231 ....: C = 2*x0*x_bar.column()
1232 ....: D = matrix.identity(x.base_ring(), n-1)
1233 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1234 ....: D = D + 2*x_bar.tensor_product(x_bar)
1235 ....: Q = matrix.block(2,2,[A,B,C,D])
1236 sage: Q == x.quadratic_representation().matrix()
1239 Test all of the properties from Theorem 11.2 in Alizadeh::
1241 sage: set_random_seed()
1242 sage: J = random_eja()
1243 sage: x,y = J.random_elements(2)
1244 sage: Lx = x.operator()
1245 sage: Lxx = (x*x).operator()
1246 sage: Qx = x.quadratic_representation()
1247 sage: Qy = y.quadratic_representation()
1248 sage: Qxy = x.quadratic_representation(y)
1249 sage: Qex = J.one().quadratic_representation(x)
1250 sage: n = ZZ.random_element(10)
1251 sage: Qxn = (x^n).quadratic_representation()
1255 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1258 Property 2 (multiply on the right for :trac:`28272`):
1260 sage: alpha = J.base_ring().random_element()
1261 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1266 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1269 sage: not x.is_invertible() or (
1272 ....: x.inverse().quadratic_representation() )
1275 sage: Qxy(J.one()) == x*y
1280 sage: not x.is_invertible() or (
1281 ....: x.quadratic_representation(x.inverse())*Qx
1282 ....: == Qx*x.quadratic_representation(x.inverse()) )
1285 sage: not x.is_invertible() or (
1286 ....: x.quadratic_representation(x.inverse())*Qx
1288 ....: 2*Lx*Qex - Qx )
1291 sage: 2*Lx*Qex - Qx == Lxx
1296 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1306 sage: not x.is_invertible() or (
1307 ....: Qx*x.inverse().operator() == Lx )
1312 sage: not x.operator_commutes_with(y) or (
1313 ....: Qx(y)^n == Qxn(y^n) )
1319 elif not other
in self
.parent():
1320 raise TypeError("'other' must live in the same algebra")
1323 M
= other
.operator()
1324 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1328 def spectral_decomposition(self
):
1330 Return the unique spectral decomposition of this element.
1334 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1335 element's left-multiplication-by operator to the subalgebra it
1336 generates. We then compute the spectral decomposition of that
1337 operator, and the spectral projectors we get back must be the
1338 left-multiplication-by operators for the idempotents we
1339 seek. Thus applying them to the identity element gives us those
1342 Since the eigenvalues are required to be distinct, we take
1343 the spectral decomposition of the zero element to be zero
1344 times the identity element of the algebra (which is idempotent,
1349 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1353 The spectral decomposition of the identity is ``1`` times itself,
1354 and the spectral decomposition of zero is ``0`` times the identity::
1356 sage: J = RealSymmetricEJA(3)
1359 sage: J.one().spectral_decomposition()
1361 sage: J.zero().spectral_decomposition()
1366 sage: J = RealSymmetricEJA(4)
1367 sage: x = sum(J.gens())
1368 sage: sd = x.spectral_decomposition()
1373 sage: c0.inner_product(c1) == 0
1375 sage: c0.is_idempotent()
1377 sage: c1.is_idempotent()
1379 sage: c0 + c1 == J.one()
1381 sage: l0*c0 + l1*c1 == x
1384 The spectral decomposition should work in subalgebras, too::
1386 sage: J = RealSymmetricEJA(4)
1387 sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
1388 sage: A = 2*b5 - 2*b8
1389 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1390 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1391 sage: (f0, f1, f2) = J1.gens()
1392 sage: f0.spectral_decomposition()
1396 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1398 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1399 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1402 def subalgebra_generated_by(self
, **kwargs
):
1404 Return the associative subalgebra of the parent EJA generated
1407 Since our parent algebra is unital, we want "subalgebra" to mean
1408 "unital subalgebra" as well; thus the subalgebra that an element
1409 generates will itself be a Euclidean Jordan algebra after
1410 restricting the algebra operations appropriately. This is the
1411 subalgebra that Faraut and Korányi work with in section II.2, for
1416 sage: from mjo.eja.eja_algebra import (random_eja,
1418 ....: RealSymmetricEJA)
1422 We can create subalgebras of Cartesian product EJAs that are not
1423 themselves Cartesian product EJAs (they're just "regular" EJAs)::
1425 sage: J1 = HadamardEJA(3)
1426 sage: J2 = RealSymmetricEJA(2)
1427 sage: J = cartesian_product([J1,J2])
1428 sage: J.one().subalgebra_generated_by()
1429 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
1433 This subalgebra, being composed of only powers, is associative::
1435 sage: set_random_seed()
1436 sage: x0 = random_eja().random_element()
1437 sage: A = x0.subalgebra_generated_by()
1438 sage: x,y,z = A.random_elements(3)
1439 sage: (x*y)*z == x*(y*z)
1442 Squaring in the subalgebra should work the same as in
1445 sage: set_random_seed()
1446 sage: x = random_eja().random_element()
1447 sage: A = x.subalgebra_generated_by()
1448 sage: A(x^2) == A(x)*A(x)
1451 By definition, the subalgebra generated by the zero element is
1452 the one-dimensional algebra generated by the identity
1453 element... unless the original algebra was trivial, in which
1454 case the subalgebra is trivial too::
1456 sage: set_random_seed()
1457 sage: A = random_eja().zero().subalgebra_generated_by()
1458 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1462 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1463 A
= self
.parent().subalgebra(powers
,
1468 A
.one
.set_cache(A(self
.parent().one()))
1472 def subalgebra_idempotent(self
):
1474 Find an idempotent in the associative subalgebra I generate
1475 using Proposition 2.3.5 in Baes.
1479 sage: from mjo.eja.eja_algebra import random_eja
1483 Ensure that we can find an idempotent in a non-trivial algebra
1484 where there are non-nilpotent elements, or that we get the dumb
1485 solution in the trivial algebra::
1487 sage: set_random_seed()
1488 sage: J = random_eja()
1489 sage: x = J.random_element()
1490 sage: while x.is_nilpotent() and not J.is_trivial():
1491 ....: x = J.random_element()
1492 sage: c = x.subalgebra_idempotent()
1497 if self
.parent().is_trivial():
1500 if self
.is_nilpotent():
1501 raise ValueError("this only works with non-nilpotent elements!")
1503 J
= self
.subalgebra_generated_by()
1506 # The image of the matrix of left-u^m-multiplication
1507 # will be minimal for some natural number s...
1509 minimal_dim
= J
.dimension()
1510 for i
in range(1, minimal_dim
):
1511 this_dim
= (u
**i
).operator().matrix().image().dimension()
1512 if this_dim
< minimal_dim
:
1513 minimal_dim
= this_dim
1516 # Now minimal_matrix should correspond to the smallest
1517 # non-zero subspace in Baes's (or really, Koecher's)
1520 # However, we need to restrict the matrix to work on the
1521 # subspace... or do we? Can't we just solve, knowing that
1522 # A(c) = u^(s+1) should have a solution in the big space,
1525 # Beware, solve_right() means that we're using COLUMN vectors.
1526 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1528 A
= u_next
.operator().matrix()
1529 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1531 # Now c is the idempotent we want, but it still lives in the subalgebra.
1532 return c
.superalgebra_element()
1537 Return my trace, the sum of my eigenvalues.
1539 In a trivial algebra, however you want to look at it, the trace is
1540 an empty sum for which we declare the result to be zero.
1544 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1551 sage: J = TrivialEJA()
1552 sage: J.zero().trace()
1556 sage: J = JordanSpinEJA(3)
1557 sage: x = sum(J.gens())
1563 sage: J = HadamardEJA(5)
1564 sage: J.one().trace()
1569 The trace of an element is a real number::
1571 sage: set_random_seed()
1572 sage: J = random_eja()
1573 sage: J.random_element().trace() in RLF
1576 The trace is linear::
1578 sage: set_random_seed()
1579 sage: J = random_eja()
1580 sage: x,y = J.random_elements(2)
1581 sage: alpha = J.base_ring().random_element()
1582 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1590 # Special case for the trivial algebra where
1591 # the trace is an empty sum.
1592 return P
.base_ring().zero()
1594 p
= P
._charpoly
_coefficients
()[r
-1]
1595 # The _charpoly_coeff function already adds the factor of
1596 # -1 to ensure that _charpoly_coeff(r-1) is really what
1597 # appears in front of t^{r-1} in the charpoly. However,
1598 # we want the negative of THAT for the trace.
1599 return -p(*self
.to_vector())
1602 def trace_inner_product(self
, other
):
1604 Return the trace inner product of myself and ``other``.
1608 sage: from mjo.eja.eja_algebra import random_eja
1612 The trace inner product is commutative, bilinear, and associative::
1614 sage: set_random_seed()
1615 sage: J = random_eja()
1616 sage: x,y,z = J.random_elements(3)
1618 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1621 sage: a = J.base_ring().random_element();
1622 sage: actual = (a*(x+z)).trace_inner_product(y)
1623 sage: expected = ( a*x.trace_inner_product(y) +
1624 ....: a*z.trace_inner_product(y) )
1625 sage: actual == expected
1627 sage: actual = x.trace_inner_product(a*(y+z))
1628 sage: expected = ( a*x.trace_inner_product(y) +
1629 ....: a*x.trace_inner_product(z) )
1630 sage: actual == expected
1633 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1637 if not other
in self
.parent():
1638 raise TypeError("'other' must live in the same algebra")
1640 return (self
*other
).trace()
1643 def trace_norm(self
):
1645 The norm of this element with respect to :meth:`trace_inner_product`.
1649 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1654 sage: J = HadamardEJA(2)
1655 sage: x = sum(J.gens())
1656 sage: x.trace_norm()
1661 sage: J = JordanSpinEJA(4)
1662 sage: x = sum(J.gens())
1663 sage: x.trace_norm()
1667 return self
.trace_inner_product(self
).sqrt()