1 from sage
.matrix
.constructor
import matrix
2 from sage
.modules
.free_module
import VectorSpace
3 from sage
.modules
.with_basis
.indexed_element
import IndexedFreeModuleElement
5 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
6 from mjo
.eja
.eja_utils
import _mat2vec
8 class FiniteDimensionalEJAElement(IndexedFreeModuleElement
):
10 An element of a Euclidean Jordan algebra.
15 Oh man, I should not be doing this. This hides the "disabled"
16 methods ``left_matrix`` and ``matrix`` from introspection;
17 in particular it removes them from tab-completion.
19 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
27 Return ``self`` raised to the power ``n``.
29 Jordan algebras are always power-associative; see for
30 example Faraut and Korányi, Proposition II.1.2 (ii).
32 We have to override this because our superclass uses row
33 vectors instead of column vectors! We, on the other hand,
34 assume column vectors everywhere.
38 sage: from mjo.eja.eja_algebra import random_eja
42 The definition of `x^2` is the unambiguous `x*x`::
44 sage: set_random_seed()
45 sage: x = random_eja().random_element()
49 A few examples of power-associativity::
51 sage: set_random_seed()
52 sage: x = random_eja().random_element()
53 sage: x*(x*x)*(x*x) == x^5
55 sage: (x*x)*(x*x*x) == x^5
58 We also know that powers operator-commute (Koecher, Chapter
61 sage: set_random_seed()
62 sage: x = random_eja().random_element()
63 sage: m = ZZ.random_element(0,10)
64 sage: n = ZZ.random_element(0,10)
65 sage: Lxm = (x^m).operator()
66 sage: Lxn = (x^n).operator()
67 sage: Lxm*Lxn == Lxn*Lxm
72 return self
.parent().one()
76 return (self
**(n
-1))*self
79 def apply_univariate_polynomial(self
, p
):
81 Apply the univariate polynomial ``p`` to this element.
83 A priori, SageMath won't allow us to apply a univariate
84 polynomial to an element of an EJA, because we don't know
85 that EJAs are rings (they are usually not associative). Of
86 course, we know that EJAs are power-associative, so the
87 operation is ultimately kosher. This function sidesteps
88 the CAS to get the answer we want and expect.
92 sage: from mjo.eja.eja_algebra import (HadamardEJA,
97 sage: R = PolynomialRing(QQ, 't')
99 sage: p = t^4 - t^3 + 5*t - 2
100 sage: J = HadamardEJA(5)
101 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
106 We should always get back an element of the algebra::
108 sage: set_random_seed()
109 sage: p = PolynomialRing(AA, 't').random_element()
110 sage: J = random_eja()
111 sage: x = J.random_element()
112 sage: x.apply_univariate_polynomial(p) in J
116 if len(p
.variables()) > 1:
117 raise ValueError("not a univariate polynomial")
120 # Convert the coeficcients to the parent's base ring,
121 # because a priori they might live in an (unnecessarily)
122 # larger ring for which P.sum() would fail below.
123 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
124 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
127 def characteristic_polynomial(self
):
129 Return the characteristic polynomial of this element.
133 sage: from mjo.eja.eja_algebra import HadamardEJA
137 The rank of `R^3` is three, and the minimal polynomial of
138 the identity element is `(t-1)` from which it follows that
139 the characteristic polynomial should be `(t-1)^3`::
141 sage: J = HadamardEJA(3)
142 sage: J.one().characteristic_polynomial()
143 t^3 - 3*t^2 + 3*t - 1
145 Likewise, the characteristic of the zero element in the
146 rank-three algebra `R^{n}` should be `t^{3}`::
148 sage: J = HadamardEJA(3)
149 sage: J.zero().characteristic_polynomial()
154 The characteristic polynomial of an element should evaluate
155 to zero on that element::
157 sage: set_random_seed()
158 sage: x = HadamardEJA(3).random_element()
159 sage: p = x.characteristic_polynomial()
160 sage: x.apply_univariate_polynomial(p)
163 The characteristic polynomials of the zero and unit elements
164 should be what we think they are in a subalgebra, too::
166 sage: J = HadamardEJA(3)
167 sage: p1 = J.one().characteristic_polynomial()
168 sage: q1 = J.zero().characteristic_polynomial()
169 sage: e0,e1,e2 = J.gens()
170 sage: A = (e0 + 2*e1 + 3*e2).subalgebra_generated_by() # dim 3
171 sage: p2 = A.one().characteristic_polynomial()
172 sage: q2 = A.zero().characteristic_polynomial()
179 p
= self
.parent().characteristic_polynomial_of()
180 return p(*self
.to_vector())
183 def inner_product(self
, other
):
185 Return the parent algebra's inner product of myself and ``other``.
189 sage: from mjo.eja.eja_algebra import (
190 ....: ComplexHermitianEJA,
192 ....: QuaternionHermitianEJA,
193 ....: RealSymmetricEJA,
198 The inner product in the Jordan spin algebra is the usual
199 inner product on `R^n` (this example only works because the
200 basis for the Jordan algebra is the standard basis in `R^n`)::
202 sage: J = JordanSpinEJA(3)
203 sage: x = vector(QQ,[1,2,3])
204 sage: y = vector(QQ,[4,5,6])
205 sage: x.inner_product(y)
207 sage: J.from_vector(x).inner_product(J.from_vector(y))
210 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
211 multiplication is the usual matrix multiplication in `S^n`,
212 so the inner product of the identity matrix with itself
215 sage: J = RealSymmetricEJA(3)
216 sage: J.one().inner_product(J.one())
219 Likewise, the inner product on `C^n` is `<X,Y> =
220 Re(trace(X*Y))`, where we must necessarily take the real
221 part because the product of Hermitian matrices may not be
224 sage: J = ComplexHermitianEJA(3)
225 sage: J.one().inner_product(J.one())
228 Ditto for the quaternions::
230 sage: J = QuaternionHermitianEJA(2)
231 sage: J.one().inner_product(J.one())
236 Ensure that we can always compute an inner product, and that
237 it gives us back a real number::
239 sage: set_random_seed()
240 sage: J = random_eja()
241 sage: x,y = J.random_elements(2)
242 sage: x.inner_product(y) in RLF
248 raise TypeError("'other' must live in the same algebra")
250 return P
.inner_product(self
, other
)
253 def operator_commutes_with(self
, other
):
255 Return whether or not this element operator-commutes
260 sage: from mjo.eja.eja_algebra import random_eja
264 The definition of a Jordan algebra says that any element
265 operator-commutes with its square::
267 sage: set_random_seed()
268 sage: x = random_eja().random_element()
269 sage: x.operator_commutes_with(x^2)
274 Test Lemma 1 from Chapter III of Koecher::
276 sage: set_random_seed()
277 sage: u,v = random_eja().random_elements(2)
278 sage: lhs = u.operator_commutes_with(u*v)
279 sage: rhs = v.operator_commutes_with(u^2)
283 Test the first polarization identity from my notes, Koecher
284 Chapter III, or from Baes (2.3)::
286 sage: set_random_seed()
287 sage: x,y = random_eja().random_elements(2)
288 sage: Lx = x.operator()
289 sage: Ly = y.operator()
290 sage: Lxx = (x*x).operator()
291 sage: Lxy = (x*y).operator()
292 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
295 Test the second polarization identity from my notes or from
298 sage: set_random_seed()
299 sage: x,y,z = random_eja().random_elements(3)
300 sage: Lx = x.operator()
301 sage: Ly = y.operator()
302 sage: Lz = z.operator()
303 sage: Lzy = (z*y).operator()
304 sage: Lxy = (x*y).operator()
305 sage: Lxz = (x*z).operator()
306 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
309 Test the third polarization identity from my notes or from
312 sage: set_random_seed()
313 sage: u,y,z = random_eja().random_elements(3)
314 sage: Lu = u.operator()
315 sage: Ly = y.operator()
316 sage: Lz = z.operator()
317 sage: Lzy = (z*y).operator()
318 sage: Luy = (u*y).operator()
319 sage: Luz = (u*z).operator()
320 sage: Luyz = (u*(y*z)).operator()
321 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
322 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
323 sage: bool(lhs == rhs)
327 if not other
in self
.parent():
328 raise TypeError("'other' must live in the same algebra")
337 Return my determinant, the product of my eigenvalues.
341 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
343 ....: RealSymmetricEJA,
344 ....: ComplexHermitianEJA,
349 sage: J = JordanSpinEJA(2)
350 sage: e0,e1 = J.gens()
351 sage: x = sum( J.gens() )
357 sage: J = JordanSpinEJA(3)
358 sage: e0,e1,e2 = J.gens()
359 sage: x = sum( J.gens() )
363 The determinant of the sole element in the rank-zero trivial
364 algebra is ``1``, by three paths of reasoning. First, its
365 characteristic polynomial is a constant ``1``, so the constant
366 term in that polynomial is ``1``. Second, the characteristic
367 polynomial evaluated at zero is again ``1``. And finally, the
368 (empty) product of its eigenvalues is likewise just unity::
370 sage: J = TrivialEJA()
376 An element is invertible if and only if its determinant is
379 sage: set_random_seed()
380 sage: x = random_eja().random_element()
381 sage: x.is_invertible() == (x.det() != 0)
384 Ensure that the determinant is multiplicative on an associative
385 subalgebra as in Faraut and Korányi's Proposition II.2.2::
387 sage: set_random_seed()
388 sage: J = random_eja().random_element().subalgebra_generated_by()
389 sage: x,y = J.random_elements(2)
390 sage: (x*y).det() == x.det()*y.det()
393 The determinant in matrix algebras is just the usual determinant::
395 sage: set_random_seed()
396 sage: X = matrix.random(QQ,3)
398 sage: J1 = RealSymmetricEJA(3)
399 sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
400 sage: expected = X.det()
401 sage: actual1 = J1(X).det()
402 sage: actual2 = J2(X).det()
403 sage: actual1 == expected
405 sage: actual2 == expected
410 sage: set_random_seed()
411 sage: J1 = ComplexHermitianEJA(2)
412 sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
413 sage: X = matrix.random(GaussianIntegers(), 2)
415 sage: expected = AA(X.det())
416 sage: actual1 = J1(J1.real_embed(X)).det()
417 sage: actual2 = J2(J2.real_embed(X)).det()
418 sage: expected == actual1
420 sage: expected == actual2
428 # Special case, since we don't get the a0=1
429 # coefficient when the rank of the algebra
431 return P
.base_ring().one()
433 p
= P
._charpoly
_coefficients
()[0]
434 # The _charpoly_coeff function already adds the factor of -1
435 # to ensure that _charpoly_coefficients()[0] is really what
436 # appears in front of t^{0} in the charpoly. However, we want
437 # (-1)^r times THAT for the determinant.
438 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 ValueError: 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 if not self
.is_invertible():
538 raise ValueError("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 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()
548 return (~self
.quadratic_representation())(self
)
551 def is_invertible(self
):
553 Return whether or not this element is invertible.
557 The usual way to do this is to check if the determinant is
558 zero, but we need the characteristic polynomial for the
559 determinant. The minimal polynomial is a lot easier to get,
560 so we use Corollary 2 in Chapter V of Koecher to check
561 whether or not the parent algebra's zero element is a root
562 of this element's minimal polynomial.
564 That is... unless the coefficients of our algebra's
565 "characteristic polynomial of" function are already cached!
566 In that case, we just use the determinant (which will be fast
569 Beware that we can't use the superclass method, because it
570 relies on the algebra being associative.
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
606 if self
.parent().is_trivial():
611 if self
.parent()._charpoly
_coefficients
.is_in_cache():
612 # The determinant will be quicker than computing the minimal
613 # polynomial from scratch, most likely.
614 return (not self
.det().is_zero())
616 # In fact, we only need to know if the constant term is non-zero,
617 # so we can pass in the field's zero element instead.
618 zero
= self
.base_ring().zero()
619 p
= self
.minimal_polynomial()
620 return not (p(zero
) == zero
)
623 def is_primitive_idempotent(self
):
625 Return whether or not this element is a primitive (or minimal)
628 A primitive idempotent is a non-zero idempotent that is not
629 the sum of two other non-zero idempotents. Remark 2.7.15 in
630 Baes shows that this is what he refers to as a "minimal
633 An element of a Euclidean Jordan algebra is a minimal idempotent
634 if it :meth:`is_idempotent` and if its Peirce subalgebra
635 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
640 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
641 ....: RealSymmetricEJA,
647 This method is sloooooow.
651 The spectral decomposition of a non-regular element should always
652 contain at least one non-minimal idempotent::
654 sage: J = RealSymmetricEJA(3)
655 sage: x = sum(J.gens())
658 sage: [ c.is_primitive_idempotent()
659 ....: for (l,c) in x.spectral_decomposition() ]
662 On the other hand, the spectral decomposition of a regular
663 element should always be in terms of minimal idempotents::
665 sage: J = JordanSpinEJA(4)
666 sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) )
669 sage: [ c.is_primitive_idempotent()
670 ....: for (l,c) in x.spectral_decomposition() ]
675 The identity element is minimal only in an EJA of rank one::
677 sage: set_random_seed()
678 sage: J = random_eja()
679 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
682 A non-idempotent cannot be a minimal idempotent::
684 sage: set_random_seed()
685 sage: J = JordanSpinEJA(4)
686 sage: x = J.random_element()
687 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
690 Proposition 2.7.19 in Baes says that an element is a minimal
691 idempotent if and only if it's idempotent with trace equal to
694 sage: set_random_seed()
695 sage: J = JordanSpinEJA(4)
696 sage: x = J.random_element()
697 sage: expected = (x.is_idempotent() and x.trace() == 1)
698 sage: actual = x.is_primitive_idempotent()
699 sage: actual == expected
702 Primitive idempotents must be non-zero::
704 sage: set_random_seed()
705 sage: J = random_eja()
706 sage: J.zero().is_idempotent()
708 sage: J.zero().is_primitive_idempotent()
711 As a consequence of the fact that primitive idempotents must
712 be non-zero, there are no primitive idempotents in a trivial
713 Euclidean Jordan algebra::
715 sage: J = TrivialEJA()
716 sage: J.one().is_idempotent()
718 sage: J.one().is_primitive_idempotent()
722 if not self
.is_idempotent():
728 (_
,_
,J1
) = self
.parent().peirce_decomposition(self
)
729 return (J1
.dimension() == 1)
732 def is_nilpotent(self
):
734 Return whether or not some power of this element is zero.
738 We use Theorem 5 in Chapter III of Koecher, which says that
739 an element ``x`` is nilpotent if and only if ``x.operator()``
740 is nilpotent. And it is a basic fact of linear algebra that
741 an operator on an `n`-dimensional space is nilpotent if and
742 only if, when raised to the `n`th power, it equals the zero
743 operator (for example, see Axler Corollary 8.8).
747 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
752 sage: J = JordanSpinEJA(3)
753 sage: x = sum(J.gens())
754 sage: x.is_nilpotent()
759 The identity element is never nilpotent, except in a trivial EJA::
761 sage: set_random_seed()
762 sage: J = random_eja()
763 sage: J.one().is_nilpotent() and not J.is_trivial()
766 The additive identity is always nilpotent::
768 sage: set_random_seed()
769 sage: random_eja().zero().is_nilpotent()
774 zero_operator
= P
.zero().operator()
775 return self
.operator()**P
.dimension() == zero_operator
778 def is_regular(self
):
780 Return whether or not this is a regular element.
784 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
789 The identity element always has degree one, but any element
790 linearly-independent from it is regular::
792 sage: J = JordanSpinEJA(5)
793 sage: J.one().is_regular()
795 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
796 sage: for x in J.gens():
797 ....: (J.one() + x).is_regular()
806 The zero element should never be regular, unless the parent
807 algebra has dimension less than or equal to one::
809 sage: set_random_seed()
810 sage: J = random_eja()
811 sage: J.dimension() <= 1 or not J.zero().is_regular()
814 The unit element isn't regular unless the algebra happens to
815 consist of only its scalar multiples::
817 sage: set_random_seed()
818 sage: J = random_eja()
819 sage: J.dimension() <= 1 or not J.one().is_regular()
823 return self
.degree() == self
.parent().rank()
828 Return the degree of this element, which is defined to be
829 the degree of its minimal polynomial.
837 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
842 sage: J = JordanSpinEJA(4)
843 sage: J.one().degree()
845 sage: e0,e1,e2,e3 = J.gens()
846 sage: (e0 - e1).degree()
849 In the spin factor algebra (of rank two), all elements that
850 aren't multiples of the identity are regular::
852 sage: set_random_seed()
853 sage: J = JordanSpinEJA.random_instance()
854 sage: n = J.dimension()
855 sage: x = J.random_element()
856 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
861 The zero and unit elements are both of degree one in nontrivial
864 sage: set_random_seed()
865 sage: J = random_eja()
866 sage: d = J.zero().degree()
867 sage: (J.is_trivial() and d == 0) or d == 1
869 sage: d = J.one().degree()
870 sage: (J.is_trivial() and d == 0) or d == 1
873 Our implementation agrees with the definition::
875 sage: set_random_seed()
876 sage: x = random_eja().random_element()
877 sage: x.degree() == x.minimal_polynomial().degree()
881 n
= self
.parent().dimension()
884 # The minimal polynomial is an empty product, i.e. the
885 # constant polynomial "1" having degree zero.
888 # The minimal polynomial of zero in a nontrivial algebra
889 # is "t", and is of degree one.
892 # If this is a nonzero element of a nontrivial algebra, it
893 # has degree at least one. It follows that, in an algebra
894 # of dimension one, the degree must be actually one.
897 # BEWARE: The subalgebra_generated_by() method uses the result
898 # of this method to construct a basis for the subalgebra. That
899 # means, in particular, that we cannot implement this method
900 # as ``self.subalgebra_generated_by().dimension()``.
902 # Algorithm: keep appending (vector representations of) powers
903 # self as rows to a matrix and echelonizing it. When its rank
904 # stops increasing, we've reached a redundancy.
906 # Given the special cases above, we can assume that "self" is
907 # nonzero, the algebra is nontrivial, and that its dimension
909 M
= matrix([(self
.parent().one()).to_vector()])
912 # Specifying the row-reduction algorithm can e.g. help over
913 # AA because it avoids the RecursionError that gets thrown
914 # when we have to look too hard for a root.
916 # Beware: QQ supports an entirely different set of "algorithm"
917 # keywords than do AA and RR.
919 from sage
.rings
.all
import QQ
920 if self
.parent().base_ring() is not QQ
:
921 algo
= "scaled_partial_pivoting"
924 M
= matrix(M
.rows() + [(self
**d
).to_vector()])
927 if new_rank
== old_rank
:
936 def left_matrix(self
):
938 Our parent class defines ``left_matrix`` and ``matrix``
939 methods whose names are misleading. We don't want them.
941 raise NotImplementedError("use operator().matrix() instead")
946 def minimal_polynomial(self
):
948 Return the minimal polynomial of this element,
949 as a function of the variable `t`.
953 We restrict ourselves to the associative subalgebra
954 generated by this element, and then return the minimal
955 polynomial of this element's operator matrix (in that
956 subalgebra). This works by Baes Proposition 2.3.16.
960 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
961 ....: RealSymmetricEJA,
967 Keeping in mind that the polynomial ``1`` evaluates the identity
968 element (also the zero element) of the trivial algebra, it is clear
969 that the polynomial ``1`` is the minimal polynomial of the only
970 element in a trivial algebra::
972 sage: J = TrivialEJA()
973 sage: J.one().minimal_polynomial()
975 sage: J.zero().minimal_polynomial()
980 The minimal polynomial of the identity and zero elements are
981 always the same, except in trivial algebras where the minimal
982 polynomial of the unit/zero element is ``1``::
984 sage: set_random_seed()
985 sage: J = random_eja()
986 sage: mu = J.one().minimal_polynomial()
987 sage: t = mu.parent().gen()
988 sage: mu + int(J.is_trivial())*(t-2)
990 sage: mu = J.zero().minimal_polynomial()
991 sage: t = mu.parent().gen()
992 sage: mu + int(J.is_trivial())*(t-1)
995 The degree of an element is (by one definition) the degree
996 of its minimal polynomial::
998 sage: set_random_seed()
999 sage: x = random_eja().random_element()
1000 sage: x.degree() == x.minimal_polynomial().degree()
1003 The minimal polynomial and the characteristic polynomial coincide
1004 and are known (see Alizadeh, Example 11.11) for all elements of
1005 the spin factor algebra that aren't scalar multiples of the
1006 identity. We require the dimension of the algebra to be at least
1007 two here so that said elements actually exist::
1009 sage: set_random_seed()
1010 sage: n_max = max(2, JordanSpinEJA._max_random_instance_size())
1011 sage: n = ZZ.random_element(2, n_max)
1012 sage: J = JordanSpinEJA(n)
1013 sage: y = J.random_element()
1014 sage: while y == y.coefficient(0)*J.one():
1015 ....: y = J.random_element()
1016 sage: y0 = y.to_vector()[0]
1017 sage: y_bar = y.to_vector()[1:]
1018 sage: actual = y.minimal_polynomial()
1019 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1020 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1021 sage: bool(actual == expected)
1024 The minimal polynomial should always kill its element::
1026 sage: set_random_seed()
1027 sage: x = random_eja().random_element()
1028 sage: p = x.minimal_polynomial()
1029 sage: x.apply_univariate_polynomial(p)
1032 The minimal polynomial is invariant under a change of basis,
1033 and in particular, a re-scaling of the basis::
1035 sage: set_random_seed()
1036 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1037 sage: n = ZZ.random_element(1, n_max)
1038 sage: J1 = RealSymmetricEJA(n)
1039 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
1040 sage: X = random_matrix(AA,n)
1041 sage: X = X*X.transpose()
1044 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
1049 # We would generate a zero-dimensional subalgebra
1050 # where the minimal polynomial would be constant.
1051 # That might be correct, but only if *this* algebra
1053 if not self
.parent().is_trivial():
1054 # Pretty sure we know what the minimal polynomial of
1055 # the zero operator is going to be. This ensures
1056 # consistency of e.g. the polynomial variable returned
1057 # in the "normal" case without us having to think about it.
1058 return self
.operator().minimal_polynomial()
1060 A
= self
.subalgebra_generated_by(orthonormalize
=False)
1061 return A(self
).operator().minimal_polynomial()
1065 def to_matrix(self
):
1067 Return an (often more natural) representation of this element as a
1070 Every finite-dimensional Euclidean Jordan Algebra is a direct
1071 sum of five simple algebras, four of which comprise Hermitian
1072 matrices. This method returns a "natural" matrix
1073 representation of this element as either a Hermitian matrix or
1078 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1079 ....: QuaternionHermitianEJA)
1083 sage: J = ComplexHermitianEJA(3)
1086 sage: J.one().to_matrix()
1096 sage: J = QuaternionHermitianEJA(2)
1099 sage: J.one().to_matrix()
1110 B
= self
.parent().matrix_basis()
1111 W
= self
.parent().matrix_space()
1113 # This is just a manual "from_vector()", but of course
1114 # matrix spaces aren't vector spaces in sage, so they
1115 # don't have a from_vector() method.
1116 return W
.linear_combination( zip(B
, self
.to_vector()) )
1121 The norm of this element with respect to :meth:`inner_product`.
1125 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1130 sage: J = HadamardEJA(2)
1131 sage: x = sum(J.gens())
1137 sage: J = JordanSpinEJA(4)
1138 sage: x = sum(J.gens())
1143 return self
.inner_product(self
).sqrt()
1148 Return the left-multiplication-by-this-element
1149 operator on the ambient algebra.
1153 sage: from mjo.eja.eja_algebra import random_eja
1157 sage: set_random_seed()
1158 sage: J = random_eja()
1159 sage: x,y = J.random_elements(2)
1160 sage: x.operator()(y) == x*y
1162 sage: y.operator()(x) == x*y
1167 left_mult_by_self
= lambda y
: self
*y
1168 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1169 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1172 def quadratic_representation(self
, other
=None):
1174 Return the quadratic representation of this element.
1178 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1183 The explicit form in the spin factor algebra is given by
1184 Alizadeh's Example 11.12::
1186 sage: set_random_seed()
1187 sage: x = JordanSpinEJA.random_instance().random_element()
1188 sage: x_vec = x.to_vector()
1189 sage: Q = matrix.identity(x.base_ring(), 0)
1190 sage: n = x_vec.degree()
1193 ....: x_bar = x_vec[1:]
1194 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1195 ....: B = 2*x0*x_bar.row()
1196 ....: C = 2*x0*x_bar.column()
1197 ....: D = matrix.identity(x.base_ring(), n-1)
1198 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1199 ....: D = D + 2*x_bar.tensor_product(x_bar)
1200 ....: Q = matrix.block(2,2,[A,B,C,D])
1201 sage: Q == x.quadratic_representation().matrix()
1204 Test all of the properties from Theorem 11.2 in Alizadeh::
1206 sage: set_random_seed()
1207 sage: J = random_eja()
1208 sage: x,y = J.random_elements(2)
1209 sage: Lx = x.operator()
1210 sage: Lxx = (x*x).operator()
1211 sage: Qx = x.quadratic_representation()
1212 sage: Qy = y.quadratic_representation()
1213 sage: Qxy = x.quadratic_representation(y)
1214 sage: Qex = J.one().quadratic_representation(x)
1215 sage: n = ZZ.random_element(10)
1216 sage: Qxn = (x^n).quadratic_representation()
1220 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1223 Property 2 (multiply on the right for :trac:`28272`):
1225 sage: alpha = J.base_ring().random_element()
1226 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1231 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1234 sage: not x.is_invertible() or (
1237 ....: x.inverse().quadratic_representation() )
1240 sage: Qxy(J.one()) == x*y
1245 sage: not x.is_invertible() or (
1246 ....: x.quadratic_representation(x.inverse())*Qx
1247 ....: == Qx*x.quadratic_representation(x.inverse()) )
1250 sage: not x.is_invertible() or (
1251 ....: x.quadratic_representation(x.inverse())*Qx
1253 ....: 2*Lx*Qex - Qx )
1256 sage: 2*Lx*Qex - Qx == Lxx
1261 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1271 sage: not x.is_invertible() or (
1272 ....: Qx*x.inverse().operator() == Lx )
1277 sage: not x.operator_commutes_with(y) or (
1278 ....: Qx(y)^n == Qxn(y^n) )
1284 elif not other
in self
.parent():
1285 raise TypeError("'other' must live in the same algebra")
1288 M
= other
.operator()
1289 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1293 def spectral_decomposition(self
):
1295 Return the unique spectral decomposition of this element.
1299 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1300 element's left-multiplication-by operator to the subalgebra it
1301 generates. We then compute the spectral decomposition of that
1302 operator, and the spectral projectors we get back must be the
1303 left-multiplication-by operators for the idempotents we
1304 seek. Thus applying them to the identity element gives us those
1307 Since the eigenvalues are required to be distinct, we take
1308 the spectral decomposition of the zero element to be zero
1309 times the identity element of the algebra (which is idempotent,
1314 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1318 The spectral decomposition of the identity is ``1`` times itself,
1319 and the spectral decomposition of zero is ``0`` times the identity::
1321 sage: J = RealSymmetricEJA(3)
1324 sage: J.one().spectral_decomposition()
1326 sage: J.zero().spectral_decomposition()
1331 sage: J = RealSymmetricEJA(4)
1332 sage: x = sum(J.gens())
1333 sage: sd = x.spectral_decomposition()
1338 sage: c0.inner_product(c1) == 0
1340 sage: c0.is_idempotent()
1342 sage: c1.is_idempotent()
1344 sage: c0 + c1 == J.one()
1346 sage: l0*c0 + l1*c1 == x
1349 The spectral decomposition should work in subalgebras, too::
1351 sage: J = RealSymmetricEJA(4)
1352 sage: (e0, e1, e2, e3, e4, e5, e6, e7, e8, e9) = J.gens()
1353 sage: A = 2*e5 - 2*e8
1354 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1355 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1356 sage: (f0, f1, f2) = J1.gens()
1357 sage: f0.spectral_decomposition()
1361 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1363 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1364 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1367 def subalgebra_generated_by(self
, **kwargs
):
1369 Return the associative subalgebra of the parent EJA generated
1372 Since our parent algebra is unital, we want "subalgebra" to mean
1373 "unital subalgebra" as well; thus the subalgebra that an element
1374 generates will itself be a Euclidean Jordan algebra after
1375 restricting the algebra operations appropriately. This is the
1376 subalgebra that Faraut and Korányi work with in section II.2, for
1381 sage: from mjo.eja.eja_algebra import random_eja
1385 This subalgebra, being composed of only powers, is associative::
1387 sage: set_random_seed()
1388 sage: x0 = random_eja().random_element()
1389 sage: A = x0.subalgebra_generated_by()
1390 sage: x,y,z = A.random_elements(3)
1391 sage: (x*y)*z == x*(y*z)
1394 Squaring in the subalgebra should work the same as in
1397 sage: set_random_seed()
1398 sage: x = random_eja().random_element()
1399 sage: A = x.subalgebra_generated_by()
1400 sage: A(x^2) == A(x)*A(x)
1403 By definition, the subalgebra generated by the zero element is
1404 the one-dimensional algebra generated by the identity
1405 element... unless the original algebra was trivial, in which
1406 case the subalgebra is trivial too::
1408 sage: set_random_seed()
1409 sage: A = random_eja().zero().subalgebra_generated_by()
1410 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1414 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1415 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1416 A
= FiniteDimensionalEJASubalgebra(self
.parent(),
1420 A
.one
.set_cache(A(self
.parent().one()))
1424 def subalgebra_idempotent(self
):
1426 Find an idempotent in the associative subalgebra I generate
1427 using Proposition 2.3.5 in Baes.
1431 sage: from mjo.eja.eja_algebra import random_eja
1435 Ensure that we can find an idempotent in a non-trivial algebra
1436 where there are non-nilpotent elements, or that we get the dumb
1437 solution in the trivial algebra::
1439 sage: set_random_seed()
1440 sage: J = random_eja()
1441 sage: x = J.random_element()
1442 sage: while x.is_nilpotent() and not J.is_trivial():
1443 ....: x = J.random_element()
1444 sage: c = x.subalgebra_idempotent()
1449 if self
.parent().is_trivial():
1452 if self
.is_nilpotent():
1453 raise ValueError("this only works with non-nilpotent elements!")
1455 J
= self
.subalgebra_generated_by()
1458 # The image of the matrix of left-u^m-multiplication
1459 # will be minimal for some natural number s...
1461 minimal_dim
= J
.dimension()
1462 for i
in range(1, minimal_dim
):
1463 this_dim
= (u
**i
).operator().matrix().image().dimension()
1464 if this_dim
< minimal_dim
:
1465 minimal_dim
= this_dim
1468 # Now minimal_matrix should correspond to the smallest
1469 # non-zero subspace in Baes's (or really, Koecher's)
1472 # However, we need to restrict the matrix to work on the
1473 # subspace... or do we? Can't we just solve, knowing that
1474 # A(c) = u^(s+1) should have a solution in the big space,
1477 # Beware, solve_right() means that we're using COLUMN vectors.
1478 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1480 A
= u_next
.operator().matrix()
1481 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1483 # Now c is the idempotent we want, but it still lives in the subalgebra.
1484 return c
.superalgebra_element()
1489 Return my trace, the sum of my eigenvalues.
1491 In a trivial algebra, however you want to look at it, the trace is
1492 an empty sum for which we declare the result to be zero.
1496 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1503 sage: J = TrivialEJA()
1504 sage: J.zero().trace()
1508 sage: J = JordanSpinEJA(3)
1509 sage: x = sum(J.gens())
1515 sage: J = HadamardEJA(5)
1516 sage: J.one().trace()
1521 The trace of an element is a real number::
1523 sage: set_random_seed()
1524 sage: J = random_eja()
1525 sage: J.random_element().trace() in RLF
1533 # Special case for the trivial algebra where
1534 # the trace is an empty sum.
1535 return P
.base_ring().zero()
1537 p
= P
._charpoly
_coefficients
()[r
-1]
1538 # The _charpoly_coeff function already adds the factor of
1539 # -1 to ensure that _charpoly_coeff(r-1) is really what
1540 # appears in front of t^{r-1} in the charpoly. However,
1541 # we want the negative of THAT for the trace.
1542 return -p(*self
.to_vector())
1545 def trace_inner_product(self
, other
):
1547 Return the trace inner product of myself and ``other``.
1551 sage: from mjo.eja.eja_algebra import random_eja
1555 The trace inner product is commutative, bilinear, and associative::
1557 sage: set_random_seed()
1558 sage: J = random_eja()
1559 sage: x,y,z = J.random_elements(3)
1561 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1564 sage: a = J.base_ring().random_element();
1565 sage: actual = (a*(x+z)).trace_inner_product(y)
1566 sage: expected = ( a*x.trace_inner_product(y) +
1567 ....: a*z.trace_inner_product(y) )
1568 sage: actual == expected
1570 sage: actual = x.trace_inner_product(a*(y+z))
1571 sage: expected = ( a*x.trace_inner_product(y) +
1572 ....: a*x.trace_inner_product(z) )
1573 sage: actual == expected
1576 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1580 if not other
in self
.parent():
1581 raise TypeError("'other' must live in the same algebra")
1583 return (self
*other
).trace()
1586 def trace_norm(self
):
1588 The norm of this element with respect to :meth:`trace_inner_product`.
1592 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1597 sage: J = HadamardEJA(2)
1598 sage: x = sum(J.gens())
1599 sage: x.trace_norm()
1604 sage: J = JordanSpinEJA(4)
1605 sage: x = sum(J.gens())
1606 sage: x.trace_norm()
1610 return self
.trace_inner_product(self
).sqrt()