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 = random_eja().random_element()
160 sage: p = x.characteristic_polynomial()
161 sage: x.apply_univariate_polynomial(p).is_zero()
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 real matrix algebras is 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
412 # Special case, since we don't get the a0=1
413 # coefficient when the rank of the algebra
415 return P
.base_ring().one()
417 p
= P
._charpoly
_coefficients
()[0]
418 # The _charpoly_coeff function already adds the factor of -1
419 # to ensure that _charpoly_coefficients()[0] is really what
420 # appears in front of t^{0} in the charpoly. However, we want
421 # (-1)^r times THAT for the determinant.
422 return ((-1)**r
)*p(*self
.to_vector())
428 Return the Jordan-multiplicative inverse of this element.
432 In general we appeal to the quadratic representation as in
433 Koecher's Theorem 12 in Chapter III, Section 5. But if the
434 parent algebra's "characteristic polynomial of" coefficients
435 happen to be cached, then we use Proposition II.2.4 in Faraut
436 and Korányi which gives a formula for the inverse based on the
437 characteristic polynomial and the Cayley-Hamilton theorem for
438 Euclidean Jordan algebras::
442 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
448 The inverse in the spin factor algebra is given in Alizadeh's
451 sage: set_random_seed()
452 sage: J = JordanSpinEJA.random_instance()
453 sage: x = J.random_element()
454 sage: while not x.is_invertible():
455 ....: x = J.random_element()
456 sage: x_vec = x.to_vector()
458 sage: x_bar = x_vec[1:]
459 sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
460 sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
461 sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
462 sage: x.inverse() == J.from_vector(x_inverse)
465 Trying to invert a non-invertible element throws an error:
467 sage: JordanSpinEJA(3).zero().inverse()
468 Traceback (most recent call last):
470 ZeroDivisionError: element is not invertible
474 The identity element is its own inverse::
476 sage: set_random_seed()
477 sage: J = random_eja()
478 sage: J.one().inverse() == J.one()
481 If an element has an inverse, it acts like one::
483 sage: set_random_seed()
484 sage: J = random_eja()
485 sage: x = J.random_element()
486 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
489 The inverse of the inverse is what we started with::
491 sage: set_random_seed()
492 sage: J = random_eja()
493 sage: x = J.random_element()
494 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
497 Proposition II.2.3 in Faraut and Korányi says that the inverse
498 of an element is the inverse of its left-multiplication operator
499 applied to the algebra's identity, when that inverse exists::
501 sage: set_random_seed()
502 sage: J = random_eja()
503 sage: x = J.random_element()
504 sage: (not x.operator().is_invertible()) or (
505 ....: x.operator().inverse()(J.one()) == x.inverse() )
508 Check that the fast (cached) and slow algorithms give the same
511 sage: set_random_seed() # long time
512 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
513 sage: x = J.random_element() # long time
514 sage: while not x.is_invertible(): # long time
515 ....: x = J.random_element() # long time
516 sage: slow = x.inverse() # long time
517 sage: _ = J._charpoly_coefficients() # long time
518 sage: fast = x.inverse() # long time
519 sage: slow == fast # long time
522 not_invertible_msg
= "element is not invertible"
523 if self
.parent()._charpoly
_coefficients
.is_in_cache():
524 # We can invert using our charpoly if it will be fast to
525 # compute. If the coefficients are cached, our rank had
527 if self
.det().is_zero():
528 raise ZeroDivisionError(not_invertible_msg
)
529 r
= self
.parent().rank()
530 a
= self
.characteristic_polynomial().coefficients(sparse
=False)
531 return (-1)**(r
+1)*sum(a
[i
+1]*self
**i
for i
in range(r
))/self
.det()
534 inv
= (~self
.quadratic_representation())(self
)
535 self
.is_invertible
.set_cache(True)
537 except ZeroDivisionError:
538 self
.is_invertible
.set_cache(False)
539 raise ZeroDivisionError(not_invertible_msg
)
543 def is_invertible(self
):
545 Return whether or not this element is invertible.
549 If computing my determinant will be fast, we do so and compare
550 with zero (Proposition II.2.4 in Faraut and
551 Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
552 reduces the problem to the invertibility of my quadratic
557 sage: from mjo.eja.eja_algebra import random_eja
561 The identity element is always invertible::
563 sage: set_random_seed()
564 sage: J = random_eja()
565 sage: J.one().is_invertible()
568 The zero element is never invertible in a non-trivial algebra::
570 sage: set_random_seed()
571 sage: J = random_eja()
572 sage: (not J.is_trivial()) and J.zero().is_invertible()
575 Test that the fast (cached) and slow algorithms give the same
578 sage: set_random_seed() # long time
579 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
580 sage: x = J.random_element() # long time
581 sage: slow = x.is_invertible() # long time
582 sage: _ = J._charpoly_coefficients() # long time
583 sage: fast = x.is_invertible() # long time
584 sage: slow == fast # long time
588 if self
.parent().is_trivial():
593 if self
.parent()._charpoly
_coefficients
.is_in_cache():
594 # The determinant will be quicker than inverting the
595 # quadratic representation, most likely.
596 return (not self
.det().is_zero())
598 # The easiest way to determine if I'm invertible is to try.
600 inv
= (~self
.quadratic_representation())(self
)
601 self
.inverse
.set_cache(inv
)
603 except ZeroDivisionError:
607 def is_primitive_idempotent(self
):
609 Return whether or not this element is a primitive (or minimal)
612 A primitive idempotent is a non-zero idempotent that is not
613 the sum of two other non-zero idempotents. Remark 2.7.15 in
614 Baes shows that this is what he refers to as a "minimal
617 An element of a Euclidean Jordan algebra is a minimal idempotent
618 if it :meth:`is_idempotent` and if its Peirce subalgebra
619 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
624 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
625 ....: RealSymmetricEJA,
631 This method is sloooooow.
635 The spectral decomposition of a non-regular element should always
636 contain at least one non-minimal idempotent::
638 sage: J = RealSymmetricEJA(3)
639 sage: x = sum(J.gens())
642 sage: [ c.is_primitive_idempotent()
643 ....: for (l,c) in x.spectral_decomposition() ]
646 On the other hand, the spectral decomposition of a regular
647 element should always be in terms of minimal idempotents::
649 sage: J = JordanSpinEJA(4)
650 sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
653 sage: [ c.is_primitive_idempotent()
654 ....: for (l,c) in x.spectral_decomposition() ]
659 The identity element is minimal only in an EJA of rank one::
661 sage: set_random_seed()
662 sage: J = random_eja()
663 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
666 A non-idempotent cannot be a minimal idempotent::
668 sage: set_random_seed()
669 sage: J = JordanSpinEJA(4)
670 sage: x = J.random_element()
671 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
674 Proposition 2.7.19 in Baes says that an element is a minimal
675 idempotent if and only if it's idempotent with trace equal to
678 sage: set_random_seed()
679 sage: J = JordanSpinEJA(4)
680 sage: x = J.random_element()
681 sage: expected = (x.is_idempotent() and x.trace() == 1)
682 sage: actual = x.is_primitive_idempotent()
683 sage: actual == expected
686 Primitive idempotents must be non-zero::
688 sage: set_random_seed()
689 sage: J = random_eja()
690 sage: J.zero().is_idempotent()
692 sage: J.zero().is_primitive_idempotent()
695 As a consequence of the fact that primitive idempotents must
696 be non-zero, there are no primitive idempotents in a trivial
697 Euclidean Jordan algebra::
699 sage: J = TrivialEJA()
700 sage: J.one().is_idempotent()
702 sage: J.one().is_primitive_idempotent()
706 if not self
.is_idempotent():
712 (_
,_
,J1
) = self
.parent().peirce_decomposition(self
)
713 return (J1
.dimension() == 1)
716 def is_nilpotent(self
):
718 Return whether or not some power of this element is zero.
722 We use Theorem 5 in Chapter III of Koecher, which says that
723 an element ``x`` is nilpotent if and only if ``x.operator()``
724 is nilpotent. And it is a basic fact of linear algebra that
725 an operator on an `n`-dimensional space is nilpotent if and
726 only if, when raised to the `n`th power, it equals the zero
727 operator (for example, see Axler Corollary 8.8).
731 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
736 sage: J = JordanSpinEJA(3)
737 sage: x = sum(J.gens())
738 sage: x.is_nilpotent()
743 The identity element is never nilpotent, except in a trivial EJA::
745 sage: set_random_seed()
746 sage: J = random_eja()
747 sage: J.one().is_nilpotent() and not J.is_trivial()
750 The additive identity is always nilpotent::
752 sage: set_random_seed()
753 sage: random_eja().zero().is_nilpotent()
758 zero_operator
= P
.zero().operator()
759 return self
.operator()**P
.dimension() == zero_operator
762 def is_regular(self
):
764 Return whether or not this is a regular element.
768 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
773 The identity element always has degree one, but any element
774 linearly-independent from it is regular::
776 sage: J = JordanSpinEJA(5)
777 sage: J.one().is_regular()
779 sage: b0, b1, b2, b3, b4 = J.gens()
782 sage: for x in J.gens():
783 ....: (J.one() + x).is_regular()
792 The zero element should never be regular, unless the parent
793 algebra has dimension less than or equal to one::
795 sage: set_random_seed()
796 sage: J = random_eja()
797 sage: J.dimension() <= 1 or not J.zero().is_regular()
800 The unit element isn't regular unless the algebra happens to
801 consist of only its scalar multiples::
803 sage: set_random_seed()
804 sage: J = random_eja()
805 sage: J.dimension() <= 1 or not J.one().is_regular()
809 return self
.degree() == self
.parent().rank()
814 Return the degree of this element, which is defined to be
815 the degree of its minimal polynomial.
823 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
828 sage: J = JordanSpinEJA(4)
829 sage: J.one().degree()
831 sage: b0,b1,b2,b3 = J.gens()
832 sage: (b0 - b1).degree()
835 In the spin factor algebra (of rank two), all elements that
836 aren't multiples of the identity are regular::
838 sage: set_random_seed()
839 sage: J = JordanSpinEJA.random_instance()
840 sage: n = J.dimension()
841 sage: x = J.random_element()
842 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
847 The zero and unit elements are both of degree one in nontrivial
850 sage: set_random_seed()
851 sage: J = random_eja()
852 sage: d = J.zero().degree()
853 sage: (J.is_trivial() and d == 0) or d == 1
855 sage: d = J.one().degree()
856 sage: (J.is_trivial() and d == 0) or d == 1
859 Our implementation agrees with the definition::
861 sage: set_random_seed()
862 sage: x = random_eja().random_element()
863 sage: x.degree() == x.minimal_polynomial().degree()
867 n
= self
.parent().dimension()
870 # The minimal polynomial is an empty product, i.e. the
871 # constant polynomial "1" having degree zero.
874 # The minimal polynomial of zero in a nontrivial algebra
875 # is "t", and is of degree one.
878 # If this is a nonzero element of a nontrivial algebra, it
879 # has degree at least one. It follows that, in an algebra
880 # of dimension one, the degree must be actually one.
883 # BEWARE: The subalgebra_generated_by() method uses the result
884 # of this method to construct a basis for the subalgebra. That
885 # means, in particular, that we cannot implement this method
886 # as ``self.subalgebra_generated_by().dimension()``.
888 # Algorithm: keep appending (vector representations of) powers
889 # self as rows to a matrix and echelonizing it. When its rank
890 # stops increasing, we've reached a redundancy.
892 # Given the special cases above, we can assume that "self" is
893 # nonzero, the algebra is nontrivial, and that its dimension
895 M
= matrix([(self
.parent().one()).to_vector()])
898 # Specifying the row-reduction algorithm can e.g. help over
899 # AA because it avoids the RecursionError that gets thrown
900 # when we have to look too hard for a root.
902 # Beware: QQ supports an entirely different set of "algorithm"
903 # keywords than do AA and RR.
905 from sage
.rings
.all
import QQ
906 if self
.parent().base_ring() is not QQ
:
907 algo
= "scaled_partial_pivoting"
910 M
= matrix(M
.rows() + [(self
**d
).to_vector()])
913 if new_rank
== old_rank
:
922 def left_matrix(self
):
924 Our parent class defines ``left_matrix`` and ``matrix``
925 methods whose names are misleading. We don't want them.
927 raise NotImplementedError("use operator().matrix() instead")
932 def minimal_polynomial(self
):
934 Return the minimal polynomial of this element,
935 as a function of the variable `t`.
939 We restrict ourselves to the associative subalgebra
940 generated by this element, and then return the minimal
941 polynomial of this element's operator matrix (in that
942 subalgebra). This works by Baes Proposition 2.3.16.
946 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
947 ....: RealSymmetricEJA,
953 Keeping in mind that the polynomial ``1`` evaluates the identity
954 element (also the zero element) of the trivial algebra, it is clear
955 that the polynomial ``1`` is the minimal polynomial of the only
956 element in a trivial algebra::
958 sage: J = TrivialEJA()
959 sage: J.one().minimal_polynomial()
961 sage: J.zero().minimal_polynomial()
966 The minimal polynomial of the identity and zero elements are
967 always the same, except in trivial algebras where the minimal
968 polynomial of the unit/zero element is ``1``::
970 sage: set_random_seed()
971 sage: J = random_eja()
972 sage: mu = J.one().minimal_polynomial()
973 sage: t = mu.parent().gen()
974 sage: mu + int(J.is_trivial())*(t-2)
976 sage: mu = J.zero().minimal_polynomial()
977 sage: t = mu.parent().gen()
978 sage: mu + int(J.is_trivial())*(t-1)
981 The degree of an element is (by one definition) the degree
982 of its minimal polynomial::
984 sage: set_random_seed()
985 sage: x = random_eja().random_element()
986 sage: x.degree() == x.minimal_polynomial().degree()
989 The minimal polynomial and the characteristic polynomial coincide
990 and are known (see Alizadeh, Example 11.11) for all elements of
991 the spin factor algebra that aren't scalar multiples of the
992 identity. We require the dimension of the algebra to be at least
993 two here so that said elements actually exist::
995 sage: set_random_seed()
996 sage: d_max = JordanSpinEJA._max_random_instance_dimension()
997 sage: n = ZZ.random_element(2, max(2,d_max))
998 sage: J = JordanSpinEJA(n)
999 sage: y = J.random_element()
1000 sage: while y == y.coefficient(0)*J.one():
1001 ....: y = J.random_element()
1002 sage: y0 = y.to_vector()[0]
1003 sage: y_bar = y.to_vector()[1:]
1004 sage: actual = y.minimal_polynomial()
1005 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1006 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1007 sage: bool(actual == expected)
1010 The minimal polynomial should always kill its element::
1012 sage: set_random_seed()
1013 sage: x = random_eja().random_element()
1014 sage: p = x.minimal_polynomial()
1015 sage: x.apply_univariate_polynomial(p)
1018 The minimal polynomial is invariant under a change of basis,
1019 and in particular, a re-scaling of the basis::
1021 sage: set_random_seed()
1022 sage: d_max = RealSymmetricEJA._max_random_instance_dimension()
1023 sage: n = ZZ.random_element(1, d_max)
1024 sage: J1 = RealSymmetricEJA(n)
1025 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
1026 sage: X = random_matrix(AA,n)
1027 sage: X = X*X.transpose()
1030 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
1035 # Pretty sure we know what the minimal polynomial of
1036 # the zero operator is going to be. This ensures
1037 # consistency of e.g. the polynomial variable returned
1038 # in the "normal" case without us having to think about it.
1039 return self
.operator().minimal_polynomial()
1041 # If we don't orthonormalize the subalgebra's basis, then the
1042 # first two monomials in the subalgebra will be self^0 and
1043 # self^1... assuming that self^1 is not a scalar multiple of
1044 # self^0 (the unit element). We special case these to avoid
1045 # having to solve a system to coerce self into the subalgebra.
1046 A
= self
.subalgebra_generated_by(orthonormalize
=False)
1048 if A
.dimension() == 1:
1049 # Does a solve to find the scalar multiple alpha such that
1050 # alpha*unit = self. We have to do this because the basis
1051 # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
1052 unit
= self
.parent().one()
1053 alpha
= self
.to_vector() / unit
.to_vector()
1054 return (unit
.operator()*alpha
).minimal_polynomial()
1056 # If the dimension of the subalgebra is >= 2, then we just
1057 # use the second basis element.
1058 return A
.monomial(1).operator().minimal_polynomial()
1062 def to_matrix(self
):
1064 Return an (often more natural) representation of this element as a
1067 Every finite-dimensional Euclidean Jordan Algebra is a direct
1068 sum of five simple algebras, four of which comprise Hermitian
1069 matrices. This method returns a "natural" matrix
1070 representation of this element as either a Hermitian matrix or
1075 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1077 ....: QuaternionHermitianEJA,
1078 ....: RealSymmetricEJA)
1082 sage: J = ComplexHermitianEJA(3)
1085 sage: J.one().to_matrix()
1096 sage: J = QuaternionHermitianEJA(2)
1099 sage: J.one().to_matrix()
1106 This also works in Cartesian product algebras::
1108 sage: J1 = HadamardEJA(1)
1109 sage: J2 = RealSymmetricEJA(2)
1110 sage: J = cartesian_product([J1,J2])
1111 sage: x = sum(J.gens())
1112 sage: x.to_matrix()[0]
1114 sage: x.to_matrix()[1]
1115 [ 1 0.7071067811865475?]
1116 [0.7071067811865475? 1]
1119 B
= self
.parent().matrix_basis()
1120 W
= self
.parent().matrix_space()
1122 if hasattr(W
, 'cartesian_factors'):
1123 # Aaaaand linear combinations don't work in Cartesian
1124 # product spaces, even though they provide a method with
1125 # that name. This is hidden behind an "if" because the
1126 # _scale() function is slow.
1127 pairs
= zip(B
, self
.to_vector())
1128 return W
.sum( _scale(b
, alpha
) for (b
,alpha
) in pairs
)
1130 # This is just a manual "from_vector()", but of course
1131 # matrix spaces aren't vector spaces in sage, so they
1132 # don't have a from_vector() method.
1133 return W
.linear_combination( zip(B
, self
.to_vector()) )
1139 The norm of this element with respect to :meth:`inner_product`.
1143 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1148 sage: J = HadamardEJA(2)
1149 sage: x = sum(J.gens())
1155 sage: J = JordanSpinEJA(4)
1156 sage: x = sum(J.gens())
1161 return self
.inner_product(self
).sqrt()
1166 Return the left-multiplication-by-this-element
1167 operator on the ambient algebra.
1171 sage: from mjo.eja.eja_algebra import random_eja
1175 sage: set_random_seed()
1176 sage: J = random_eja()
1177 sage: x,y = J.random_elements(2)
1178 sage: x.operator()(y) == x*y
1180 sage: y.operator()(x) == x*y
1185 left_mult_by_self
= lambda y
: self
*y
1186 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1187 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1190 def quadratic_representation(self
, other
=None):
1192 Return the quadratic representation of this element.
1196 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1201 The explicit form in the spin factor algebra is given by
1202 Alizadeh's Example 11.12::
1204 sage: set_random_seed()
1205 sage: x = JordanSpinEJA.random_instance().random_element()
1206 sage: x_vec = x.to_vector()
1207 sage: Q = matrix.identity(x.base_ring(), 0)
1208 sage: n = x_vec.degree()
1211 ....: x_bar = x_vec[1:]
1212 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1213 ....: B = 2*x0*x_bar.row()
1214 ....: C = 2*x0*x_bar.column()
1215 ....: D = matrix.identity(x.base_ring(), n-1)
1216 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1217 ....: D = D + 2*x_bar.tensor_product(x_bar)
1218 ....: Q = matrix.block(2,2,[A,B,C,D])
1219 sage: Q == x.quadratic_representation().matrix()
1222 Test all of the properties from Theorem 11.2 in Alizadeh::
1224 sage: set_random_seed()
1225 sage: J = random_eja()
1226 sage: x,y = J.random_elements(2)
1227 sage: Lx = x.operator()
1228 sage: Lxx = (x*x).operator()
1229 sage: Qx = x.quadratic_representation()
1230 sage: Qy = y.quadratic_representation()
1231 sage: Qxy = x.quadratic_representation(y)
1232 sage: Qex = J.one().quadratic_representation(x)
1233 sage: n = ZZ.random_element(10)
1234 sage: Qxn = (x^n).quadratic_representation()
1238 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1241 Property 2 (multiply on the right for :trac:`28272`):
1243 sage: alpha = J.base_ring().random_element()
1244 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1249 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1252 sage: not x.is_invertible() or (
1255 ....: x.inverse().quadratic_representation() )
1258 sage: Qxy(J.one()) == x*y
1263 sage: not x.is_invertible() or (
1264 ....: x.quadratic_representation(x.inverse())*Qx
1265 ....: == Qx*x.quadratic_representation(x.inverse()) )
1268 sage: not x.is_invertible() or (
1269 ....: x.quadratic_representation(x.inverse())*Qx
1271 ....: 2*Lx*Qex - Qx )
1274 sage: 2*Lx*Qex - Qx == Lxx
1279 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1289 sage: not x.is_invertible() or (
1290 ....: Qx*x.inverse().operator() == Lx )
1295 sage: not x.operator_commutes_with(y) or (
1296 ....: Qx(y)^n == Qxn(y^n) )
1302 elif not other
in self
.parent():
1303 raise TypeError("'other' must live in the same algebra")
1306 M
= other
.operator()
1307 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1311 def spectral_decomposition(self
):
1313 Return the unique spectral decomposition of this element.
1317 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1318 element's left-multiplication-by operator to the subalgebra it
1319 generates. We then compute the spectral decomposition of that
1320 operator, and the spectral projectors we get back must be the
1321 left-multiplication-by operators for the idempotents we
1322 seek. Thus applying them to the identity element gives us those
1325 Since the eigenvalues are required to be distinct, we take
1326 the spectral decomposition of the zero element to be zero
1327 times the identity element of the algebra (which is idempotent,
1332 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1336 The spectral decomposition of the identity is ``1`` times itself,
1337 and the spectral decomposition of zero is ``0`` times the identity::
1339 sage: J = RealSymmetricEJA(3)
1342 sage: J.one().spectral_decomposition()
1344 sage: J.zero().spectral_decomposition()
1349 sage: J = RealSymmetricEJA(4)
1350 sage: x = sum(J.gens())
1351 sage: sd = x.spectral_decomposition()
1356 sage: c0.inner_product(c1) == 0
1358 sage: c0.is_idempotent()
1360 sage: c1.is_idempotent()
1362 sage: c0 + c1 == J.one()
1364 sage: l0*c0 + l1*c1 == x
1367 The spectral decomposition should work in subalgebras, too::
1369 sage: J = RealSymmetricEJA(4)
1370 sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
1371 sage: A = 2*b5 - 2*b8
1372 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1373 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1374 sage: (f0, f1, f2) = J1.gens()
1375 sage: f0.spectral_decomposition()
1379 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1381 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1382 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1385 def subalgebra_generated_by(self
, **kwargs
):
1387 Return the associative subalgebra of the parent EJA generated
1390 Since our parent algebra is unital, we want "subalgebra" to mean
1391 "unital subalgebra" as well; thus the subalgebra that an element
1392 generates will itself be a Euclidean Jordan algebra after
1393 restricting the algebra operations appropriately. This is the
1394 subalgebra that Faraut and Korányi work with in section II.2, for
1399 sage: from mjo.eja.eja_algebra import (random_eja,
1401 ....: RealSymmetricEJA)
1405 We can create subalgebras of Cartesian product EJAs that are not
1406 themselves Cartesian product EJAs (they're just "regular" EJAs)::
1408 sage: J1 = HadamardEJA(3)
1409 sage: J2 = RealSymmetricEJA(2)
1410 sage: J = cartesian_product([J1,J2])
1411 sage: J.one().subalgebra_generated_by()
1412 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
1416 This subalgebra, being composed of only powers, is associative::
1418 sage: set_random_seed()
1419 sage: x0 = random_eja().random_element()
1420 sage: A = x0.subalgebra_generated_by()
1421 sage: x,y,z = A.random_elements(3)
1422 sage: (x*y)*z == x*(y*z)
1425 Squaring in the subalgebra should work the same as in
1428 sage: set_random_seed()
1429 sage: x = random_eja().random_element()
1430 sage: A = x.subalgebra_generated_by()
1431 sage: A(x^2) == A(x)*A(x)
1434 By definition, the subalgebra generated by the zero element is
1435 the one-dimensional algebra generated by the identity
1436 element... unless the original algebra was trivial, in which
1437 case the subalgebra is trivial too::
1439 sage: set_random_seed()
1440 sage: A = random_eja().zero().subalgebra_generated_by()
1441 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1445 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1446 A
= self
.parent().subalgebra(powers
,
1451 A
.one
.set_cache(A(self
.parent().one()))
1455 def subalgebra_idempotent(self
):
1457 Find an idempotent in the associative subalgebra I generate
1458 using Proposition 2.3.5 in Baes.
1462 sage: from mjo.eja.eja_algebra import random_eja
1466 Ensure that we can find an idempotent in a non-trivial algebra
1467 where there are non-nilpotent elements, or that we get the dumb
1468 solution in the trivial algebra::
1470 sage: set_random_seed()
1471 sage: J = random_eja()
1472 sage: x = J.random_element()
1473 sage: while x.is_nilpotent() and not J.is_trivial():
1474 ....: x = J.random_element()
1475 sage: c = x.subalgebra_idempotent()
1480 if self
.parent().is_trivial():
1483 if self
.is_nilpotent():
1484 raise ValueError("this only works with non-nilpotent elements!")
1486 J
= self
.subalgebra_generated_by()
1489 # The image of the matrix of left-u^m-multiplication
1490 # will be minimal for some natural number s...
1492 minimal_dim
= J
.dimension()
1493 for i
in range(1, minimal_dim
):
1494 this_dim
= (u
**i
).operator().matrix().image().dimension()
1495 if this_dim
< minimal_dim
:
1496 minimal_dim
= this_dim
1499 # Now minimal_matrix should correspond to the smallest
1500 # non-zero subspace in Baes's (or really, Koecher's)
1503 # However, we need to restrict the matrix to work on the
1504 # subspace... or do we? Can't we just solve, knowing that
1505 # A(c) = u^(s+1) should have a solution in the big space,
1508 # Beware, solve_right() means that we're using COLUMN vectors.
1509 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1511 A
= u_next
.operator().matrix()
1512 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1514 # Now c is the idempotent we want, but it still lives in the subalgebra.
1515 return c
.superalgebra_element()
1520 Return my trace, the sum of my eigenvalues.
1522 In a trivial algebra, however you want to look at it, the trace is
1523 an empty sum for which we declare the result to be zero.
1527 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1534 sage: J = TrivialEJA()
1535 sage: J.zero().trace()
1539 sage: J = JordanSpinEJA(3)
1540 sage: x = sum(J.gens())
1546 sage: J = HadamardEJA(5)
1547 sage: J.one().trace()
1552 The trace of an element is a real number::
1554 sage: set_random_seed()
1555 sage: J = random_eja()
1556 sage: J.random_element().trace() in RLF
1559 The trace is linear::
1561 sage: set_random_seed()
1562 sage: J = random_eja()
1563 sage: x,y = J.random_elements(2)
1564 sage: alpha = J.base_ring().random_element()
1565 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1573 # Special case for the trivial algebra where
1574 # the trace is an empty sum.
1575 return P
.base_ring().zero()
1577 p
= P
._charpoly
_coefficients
()[r
-1]
1578 # The _charpoly_coeff function already adds the factor of
1579 # -1 to ensure that _charpoly_coeff(r-1) is really what
1580 # appears in front of t^{r-1} in the charpoly. However,
1581 # we want the negative of THAT for the trace.
1582 return -p(*self
.to_vector())
1585 def trace_inner_product(self
, other
):
1587 Return the trace inner product of myself and ``other``.
1591 sage: from mjo.eja.eja_algebra import random_eja
1595 The trace inner product is commutative, bilinear, and associative::
1597 sage: set_random_seed()
1598 sage: J = random_eja()
1599 sage: x,y,z = J.random_elements(3)
1601 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1604 sage: a = J.base_ring().random_element();
1605 sage: actual = (a*(x+z)).trace_inner_product(y)
1606 sage: expected = ( a*x.trace_inner_product(y) +
1607 ....: a*z.trace_inner_product(y) )
1608 sage: actual == expected
1610 sage: actual = x.trace_inner_product(a*(y+z))
1611 sage: expected = ( a*x.trace_inner_product(y) +
1612 ....: a*x.trace_inner_product(z) )
1613 sage: actual == expected
1616 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1620 if not other
in self
.parent():
1621 raise TypeError("'other' must live in the same algebra")
1623 return (self
*other
).trace()
1626 def trace_norm(self
):
1628 The norm of this element with respect to :meth:`trace_inner_product`.
1632 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1637 sage: J = HadamardEJA(2)
1638 sage: x = sum(J.gens())
1639 sage: x.trace_norm()
1644 sage: J = JordanSpinEJA(4)
1645 sage: x = sum(J.gens())
1646 sage: x.trace_norm()
1650 return self
.trace_inner_product(self
).sqrt()