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 EJAOperator
7 from mjo
.eja
.eja_utils
import _scale
10 class EJAElement(IndexedFreeModuleElement
):
12 An element of a Euclidean Jordan algebra.
17 Oh man, I should not be doing this. This hides the "disabled"
18 methods ``left_matrix`` and ``matrix`` from introspection;
19 in particular it removes them from tab-completion.
21 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
29 Return ``self`` raised to the power ``n``.
31 Jordan algebras are always power-associative; see for
32 example Faraut and Korányi, Proposition II.1.2 (ii).
34 We have to override this because our superclass uses row
35 vectors instead of column vectors! We, on the other hand,
36 assume column vectors everywhere.
40 sage: from mjo.eja.eja_algebra import random_eja
44 The definition of `x^2` is the unambiguous `x*x`::
46 sage: x = random_eja().random_element()
50 A few examples of power-associativity::
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: x = random_eja().random_element()
62 sage: m = ZZ.random_element(0,10)
63 sage: n = ZZ.random_element(0,10)
64 sage: Lxm = (x^m).operator()
65 sage: Lxn = (x^n).operator()
66 sage: Lxm*Lxn == Lxn*Lxm
71 return self
.parent().one()
75 return (self
**(n
-1))*self
78 def apply_univariate_polynomial(self
, p
):
80 Apply the univariate polynomial ``p`` to this element.
82 A priori, SageMath won't allow us to apply a univariate
83 polynomial to an element of an EJA, because we don't know
84 that EJAs are rings (they are usually not associative). Of
85 course, we know that EJAs are power-associative, so the
86 operation is ultimately kosher. This function sidesteps
87 the CAS to get the answer we want and expect.
91 sage: from mjo.eja.eja_algebra import (HadamardEJA,
96 sage: R = PolynomialRing(QQ, 't')
98 sage: p = t^4 - t^3 + 5*t - 2
99 sage: J = HadamardEJA(5)
100 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
105 We should always get back an element of the algebra::
107 sage: p = PolynomialRing(AA, 't').random_element()
108 sage: J = random_eja()
109 sage: x = J.random_element()
110 sage: x.apply_univariate_polynomial(p) in J
114 if len(p
.variables()) > 1:
115 raise ValueError("not a univariate polynomial")
118 # Convert the coeficcients to the parent's base ring,
119 # because a priori they might live in an (unnecessarily)
120 # larger ring for which P.sum() would fail below.
121 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
122 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
125 def characteristic_polynomial(self
):
127 Return the characteristic polynomial of this element.
131 sage: from mjo.eja.eja_algebra import (random_eja,
136 The rank of `R^3` is three, and the minimal polynomial of
137 the identity element is `(t-1)` from which it follows that
138 the characteristic polynomial should be `(t-1)^3`::
140 sage: J = HadamardEJA(3)
141 sage: J.one().characteristic_polynomial()
142 t^3 - 3*t^2 + 3*t - 1
144 Likewise, the characteristic of the zero element in the
145 rank-three algebra `R^{n}` should be `t^{3}`::
147 sage: J = HadamardEJA(3)
148 sage: J.zero().characteristic_polynomial()
153 The characteristic polynomial of an element should evaluate
154 to zero on that element::
156 sage: x = random_eja().random_element()
157 sage: p = x.characteristic_polynomial()
158 sage: x.apply_univariate_polynomial(p).is_zero()
161 The characteristic polynomials of the zero and unit elements
162 should be what we think they are in a subalgebra, too::
164 sage: J = HadamardEJA(3)
165 sage: p1 = J.one().characteristic_polynomial()
166 sage: q1 = J.zero().characteristic_polynomial()
167 sage: b0,b1,b2 = J.gens()
168 sage: A = (b0 + 2*b1 + 3*b2).subalgebra_generated_by() # dim 3
169 sage: p2 = A.one().characteristic_polynomial()
170 sage: q2 = A.zero().characteristic_polynomial()
177 p
= self
.parent().characteristic_polynomial_of()
178 return p(*self
.to_vector())
181 def inner_product(self
, other
):
183 Return the parent algebra's inner product of myself and ``other``.
187 sage: from mjo.eja.eja_algebra import (
188 ....: ComplexHermitianEJA,
190 ....: QuaternionHermitianEJA,
191 ....: RealSymmetricEJA,
196 The inner product in the Jordan spin algebra is the usual
197 inner product on `R^n` (this example only works because the
198 basis for the Jordan algebra is the standard basis in `R^n`)::
200 sage: J = JordanSpinEJA(3)
201 sage: x = vector(QQ,[1,2,3])
202 sage: y = vector(QQ,[4,5,6])
203 sage: x.inner_product(y)
205 sage: J.from_vector(x).inner_product(J.from_vector(y))
208 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
209 multiplication is the usual matrix multiplication in `S^n`,
210 so the inner product of the identity matrix with itself
213 sage: J = RealSymmetricEJA(3)
214 sage: J.one().inner_product(J.one())
217 Likewise, the inner product on `C^n` is `<X,Y> =
218 Re(trace(X*Y))`, where we must necessarily take the real
219 part because the product of Hermitian matrices may not be
222 sage: J = ComplexHermitianEJA(3)
223 sage: J.one().inner_product(J.one())
226 Ditto for the quaternions::
228 sage: J = QuaternionHermitianEJA(2)
229 sage: J.one().inner_product(J.one())
234 Ensure that we can always compute an inner product, and that
235 it gives us back a real number::
237 sage: J = random_eja()
238 sage: x,y = J.random_elements(2)
239 sage: x.inner_product(y) in RLF
245 raise TypeError("'other' must live in the same algebra")
247 return P
.inner_product(self
, other
)
250 def operator_commutes_with(self
, other
):
252 Return whether or not this element operator-commutes
257 sage: from mjo.eja.eja_algebra import random_eja
261 The definition of a Jordan algebra says that any element
262 operator-commutes with its square::
264 sage: x = random_eja().random_element()
265 sage: x.operator_commutes_with(x^2)
270 Test Lemma 1 from Chapter III of Koecher::
272 sage: u,v = random_eja().random_elements(2)
273 sage: lhs = u.operator_commutes_with(u*v)
274 sage: rhs = v.operator_commutes_with(u^2)
278 Test the first polarization identity from my notes, Koecher
279 Chapter III, or from Baes (2.3)::
281 sage: x,y = random_eja().random_elements(2)
282 sage: Lx = x.operator()
283 sage: Ly = y.operator()
284 sage: Lxx = (x*x).operator()
285 sage: Lxy = (x*y).operator()
286 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
289 Test the second polarization identity from my notes or from
292 sage: x,y,z = random_eja().random_elements(3) # long time
293 sage: Lx = x.operator() # long time
294 sage: Ly = y.operator() # long time
295 sage: Lz = z.operator() # long time
296 sage: Lzy = (z*y).operator() # long time
297 sage: Lxy = (x*y).operator() # long time
298 sage: Lxz = (x*z).operator() # long time
299 sage: lhs = Lx*Lzy + Lz*Lxy + Ly*Lxz # long time
300 sage: rhs = Lzy*Lx + Lxy*Lz + Lxz*Ly # long time
301 sage: bool(lhs == rhs) # long time
304 Test the third polarization identity from my notes or from
307 sage: u,y,z = random_eja().random_elements(3) # long time
308 sage: Lu = u.operator() # long time
309 sage: Ly = y.operator() # long time
310 sage: Lz = z.operator() # long time
311 sage: Lzy = (z*y).operator() # long time
312 sage: Luy = (u*y).operator() # long time
313 sage: Luz = (u*z).operator() # long time
314 sage: Luyz = (u*(y*z)).operator() # long time
315 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz # long time
316 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly # long time
317 sage: bool(lhs == rhs) # long time
321 if not other
in self
.parent():
322 raise TypeError("'other' must live in the same algebra")
331 Return my determinant, the product of my eigenvalues.
335 sage: from mjo.eja.eja_algebra import (AlbertEJA,
338 ....: RealSymmetricEJA,
339 ....: ComplexHermitianEJA,
344 sage: J = JordanSpinEJA(2)
345 sage: x = sum( J.gens() )
351 sage: J = JordanSpinEJA(3)
352 sage: x = sum( J.gens() )
356 The determinant of the sole element in the rank-zero trivial
357 algebra is ``1``, by three paths of reasoning. First, its
358 characteristic polynomial is a constant ``1``, so the constant
359 term in that polynomial is ``1``. Second, the characteristic
360 polynomial evaluated at zero is again ``1``. And finally, the
361 (empty) product of its eigenvalues is likewise just unity::
363 sage: J = TrivialEJA()
369 An element is invertible if and only if its determinant is
372 sage: x = random_eja().random_element()
373 sage: x.is_invertible() == (x.det() != 0)
376 Ensure that the determinant is multiplicative on an associative
377 subalgebra as in Faraut and Korányi's Proposition II.2.2::
379 sage: x0 = random_eja().random_element()
380 sage: J = x0.subalgebra_generated_by(orthonormalize=False)
381 sage: x,y = J.random_elements(2)
382 sage: (x*y).det() == x.det()*y.det()
385 The determinant in real matrix algebras is the usual determinant::
387 sage: X = matrix.random(QQ,3)
389 sage: J1 = RealSymmetricEJA(3)
390 sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
391 sage: expected = X.det()
392 sage: actual1 = J1(X).det()
393 sage: actual2 = J2(X).det()
394 sage: actual1 == expected
396 sage: actual2 == expected
399 There's a formula for the determinant of the Albert algebra
400 (Yokota, Section 2.1)::
402 sage: def albert_det(x):
403 ....: X = x.to_matrix()
404 ....: res = X[0,0]*X[1,1]*X[2,2]
405 ....: res += 2*(X[1,2]*X[2,0]*X[0,1]).real()
406 ....: res -= X[0,0]*X[1,2]*X[2,1]
407 ....: res -= X[1,1]*X[2,0]*X[0,2]
408 ....: res -= X[2,2]*X[0,1]*X[1,0]
409 ....: return res.leading_coefficient()
410 sage: J = AlbertEJA(field=QQ, orthonormalize=False)
411 sage: xs = J.random_elements(10)
412 sage: all( albert_det(x) == x.det() for x in xs )
420 # Special case, since we don't get the a0=1
421 # coefficient when the rank of the algebra
423 return P
.base_ring().one()
425 p
= P
._charpoly
_coefficients
()[0]
426 # The _charpoly_coeff function already adds the factor of -1
427 # to ensure that _charpoly_coefficients()[0] is really what
428 # appears in front of t^{0} in the charpoly. However, we want
429 # (-1)^r times THAT for the determinant.
430 return ((-1)**r
)*p(*self
.to_vector())
436 Return the Jordan-multiplicative inverse of this element.
440 In general we appeal to the quadratic representation as in
441 Koecher's Theorem 12 in Chapter III, Section 5. But if the
442 parent algebra's "characteristic polynomial of" coefficients
443 happen to be cached, then we use Proposition II.2.4 in Faraut
444 and Korányi which gives a formula for the inverse based on the
445 characteristic polynomial and the Cayley-Hamilton theorem for
446 Euclidean Jordan algebras::
450 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
456 The inverse in the spin factor algebra is given in Alizadeh's
459 sage: J = JordanSpinEJA.random_instance()
460 sage: x = J.random_element()
461 sage: while not x.is_invertible():
462 ....: x = J.random_element()
463 sage: x_vec = x.to_vector()
465 sage: x_bar = x_vec[1:]
466 sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
467 sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
468 sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
469 sage: x.inverse() == J.from_vector(x_inverse)
472 Trying to invert a non-invertible element throws an error:
474 sage: JordanSpinEJA(3).zero().inverse()
475 Traceback (most recent call last):
477 ZeroDivisionError: element is not invertible
481 The identity element is its own inverse::
483 sage: J = random_eja()
484 sage: J.one().inverse() == J.one()
487 If an element has an inverse, it acts like one::
489 sage: J = random_eja()
490 sage: x = J.random_element()
491 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
494 The inverse of the inverse is what we started with::
496 sage: J = random_eja()
497 sage: x = J.random_element()
498 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
501 Proposition II.2.3 in Faraut and Korányi says that the inverse
502 of an element is the inverse of its left-multiplication operator
503 applied to the algebra's identity, when that inverse exists::
505 sage: J = random_eja() # long time
506 sage: x = J.random_element() # long time
507 sage: (not x.operator().is_invertible()) or ( # long time
508 ....: x.operator().inverse()(J.one()) # long time
510 ....: x.inverse() ) # long time
513 Check that the fast (cached) and slow algorithms give the same
516 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
517 sage: x = J.random_element() # long time
518 sage: while not x.is_invertible(): # long time
519 ....: x = J.random_element() # long time
520 sage: slow = x.inverse() # long time
521 sage: _ = J._charpoly_coefficients() # long time
522 sage: fast = x.inverse() # long time
523 sage: slow == fast # long time
526 not_invertible_msg
= "element is not invertible"
528 algebra
= self
.parent()
529 if algebra
._charpoly
_coefficients
.is_in_cache():
530 # We can invert using our charpoly if it will be fast to
531 # compute. If the coefficients are cached, our rank had
533 if self
.det().is_zero():
534 raise ZeroDivisionError(not_invertible_msg
)
536 a
= self
.characteristic_polynomial().coefficients(sparse
=False)
537 return (-1)**(r
+1)*algebra
.sum(a
[i
+1]*self
**i
538 for i
in range(r
))/self
.det()
541 inv
= (~self
.quadratic_representation())(self
)
542 self
.is_invertible
.set_cache(True)
544 except ZeroDivisionError:
545 self
.is_invertible
.set_cache(False)
546 raise ZeroDivisionError(not_invertible_msg
)
550 def is_invertible(self
):
552 Return whether or not this element is invertible.
556 If computing my determinant will be fast, we do so and compare
557 with zero (Proposition II.2.4 in Faraut and
558 Korányi). Otherwise, Proposition II.3.2 in Faraut and Korányi
559 reduces the problem to the invertibility of my quadratic
564 sage: from mjo.eja.eja_algebra import random_eja
568 The identity element is always invertible::
570 sage: J = random_eja()
571 sage: J.one().is_invertible()
574 The zero element is never invertible in a non-trivial algebra::
576 sage: J = random_eja()
577 sage: (not J.is_trivial()) and J.zero().is_invertible()
580 Test that the fast (cached) and slow algorithms give the same
583 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
584 sage: x = J.random_element() # long time
585 sage: slow = x.is_invertible() # long time
586 sage: _ = J._charpoly_coefficients() # long time
587 sage: fast = x.is_invertible() # long time
588 sage: slow == fast # long time
592 if self
.parent().is_trivial():
597 if self
.parent()._charpoly
_coefficients
.is_in_cache():
598 # The determinant will be quicker than inverting the
599 # quadratic representation, most likely.
600 return (not self
.det().is_zero())
602 # The easiest way to determine if I'm invertible is to try.
604 inv
= (~self
.quadratic_representation())(self
)
605 self
.inverse
.set_cache(inv
)
607 except ZeroDivisionError:
611 def is_primitive_idempotent(self
):
613 Return whether or not this element is a primitive (or minimal)
616 A primitive idempotent is a non-zero idempotent that is not
617 the sum of two other non-zero idempotents. Remark 2.7.15 in
618 Baes shows that this is what he refers to as a "minimal
621 An element of a Euclidean Jordan algebra is a minimal idempotent
622 if it :meth:`is_idempotent` and if its Peirce subalgebra
623 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
628 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
629 ....: RealSymmetricEJA,
635 This method is sloooooow.
639 The spectral decomposition of a non-regular element should always
640 contain at least one non-minimal idempotent::
642 sage: J = RealSymmetricEJA(3)
643 sage: x = sum(J.gens())
646 sage: [ c.is_primitive_idempotent()
647 ....: for (l,c) in x.spectral_decomposition() ]
650 On the other hand, the spectral decomposition of a regular
651 element should always be in terms of minimal idempotents::
653 sage: J = JordanSpinEJA(4)
654 sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
657 sage: [ c.is_primitive_idempotent()
658 ....: for (l,c) in x.spectral_decomposition() ]
663 The identity element is minimal only in an EJA of rank one::
665 sage: J = random_eja()
666 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
669 A non-idempotent cannot be a minimal idempotent::
671 sage: J = JordanSpinEJA(4)
672 sage: x = J.random_element()
673 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
676 Proposition 2.7.19 in Baes says that an element is a minimal
677 idempotent if and only if it's idempotent with trace equal to
680 sage: J = JordanSpinEJA(4)
681 sage: x = J.random_element()
682 sage: expected = (x.is_idempotent() and x.trace() == 1)
683 sage: actual = x.is_primitive_idempotent()
684 sage: actual == expected
687 Primitive idempotents must be non-zero::
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: J = random_eja()
746 sage: J.one().is_nilpotent() and not J.is_trivial()
749 The additive identity is always nilpotent::
751 sage: random_eja().zero().is_nilpotent()
756 zero_operator
= P
.zero().operator()
757 return self
.operator()**P
.dimension() == zero_operator
760 def is_regular(self
):
762 Return whether or not this is a regular element.
766 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
771 The identity element always has degree one, but any element
772 linearly-independent from it is regular::
774 sage: J = JordanSpinEJA(5)
775 sage: J.one().is_regular()
777 sage: b0, b1, b2, b3, b4 = J.gens()
780 sage: for x in J.gens():
781 ....: (J.one() + x).is_regular()
790 The zero element should never be regular, unless the parent
791 algebra has dimension less than or equal to one::
793 sage: J = random_eja()
794 sage: J.dimension() <= 1 or not J.zero().is_regular()
797 The unit element isn't regular unless the algebra happens to
798 consist of only its scalar multiples::
800 sage: J = random_eja()
801 sage: J.dimension() <= 1 or not J.one().is_regular()
805 return self
.degree() == self
.parent().rank()
810 Return the degree of this element, which is defined to be
811 the degree of its minimal polynomial.
815 First we handle the special cases where the algebra is
816 trivial, this element is zero, or the dimension of the algebra
817 is one and this element is not zero. With those out of the
818 way, we may assume that ``self`` is nonzero, the algebra is
819 nontrivial, and that the dimension of the algebra is at least
822 Beginning with the algebra's unit element (power zero), we add
823 successive (basis representations of) powers of this element
824 to a matrix, row-reducing at each step. After row-reducing, we
825 check the rank of the matrix. If adding a row and row-reducing
826 does not increase the rank of the matrix at any point, the row
827 we've just added lives in the span of the previous ones; thus
828 the corresponding power of ``self`` lives in the span of its
829 lesser powers. When that happens, the degree of the minimal
830 polynomial is the rank of the matrix; if it never happens, the
831 degree must be the dimension of the entire space.
835 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
840 sage: J = JordanSpinEJA(4)
841 sage: J.one().degree()
843 sage: b0,b1,b2,b3 = J.gens()
844 sage: (b0 - b1).degree()
847 In the spin factor algebra (of rank two), all elements that
848 aren't multiples of the identity are regular::
850 sage: J = JordanSpinEJA.random_instance()
851 sage: n = J.dimension()
852 sage: x = J.random_element()
853 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
858 The zero and unit elements are both of degree one in nontrivial
861 sage: J = random_eja()
862 sage: d = J.zero().degree()
863 sage: (J.is_trivial() and d == 0) or d == 1
865 sage: d = J.one().degree()
866 sage: (J.is_trivial() and d == 0) or d == 1
869 Our implementation agrees with the definition::
871 sage: x = random_eja().random_element()
872 sage: x.degree() == x.minimal_polynomial().degree()
875 n
= self
.parent().dimension()
878 # The minimal polynomial is an empty product, i.e. the
879 # constant polynomial "1" having degree zero.
882 # The minimal polynomial of zero in a nontrivial algebra
883 # is "t", and is of degree one.
886 # If this is a nonzero element of a nontrivial algebra, it
887 # has degree at least one. It follows that, in an algebra
888 # of dimension one, the degree must be actually one.
891 # BEWARE: The subalgebra_generated_by() method uses the result
892 # of this method to construct a basis for the subalgebra. That
893 # means, in particular, that we cannot implement this method
894 # as ``self.subalgebra_generated_by().dimension()``.
896 # Algorithm: keep appending (vector representations of) powers
897 # self as rows to a matrix and echelonizing it. When its rank
898 # stops increasing, we've reached a redundancy.
900 # Given the special cases above, we can assume that "self" is
901 # nonzero, the algebra is nontrivial, and that its dimension
903 M
= matrix([(self
.parent().one()).to_vector()])
906 # Specifying the row-reduction algorithm can e.g. help over
907 # AA because it avoids the RecursionError that gets thrown
908 # when we have to look too hard for a root.
910 # Beware: QQ supports an entirely different set of "algorithm"
911 # keywords than do AA and RR.
913 from sage
.rings
.all
import QQ
914 if self
.parent().base_ring() is not QQ
:
915 algo
= "scaled_partial_pivoting"
918 M
= matrix(M
.rows() + [(self
**d
).to_vector()])
921 if new_rank
== old_rank
:
930 def left_matrix(self
):
932 Our parent class defines ``left_matrix`` and ``matrix``
933 methods whose names are misleading. We don't want them.
935 raise NotImplementedError("use operator().matrix() instead")
940 def minimal_polynomial(self
):
942 Return the minimal polynomial of this element,
943 as a function of the variable `t`.
947 We restrict ourselves to the associative subalgebra
948 generated by this element, and then return the minimal
949 polynomial of this element's operator matrix (in that
950 subalgebra). This works by Baes Proposition 2.3.16.
954 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
955 ....: RealSymmetricEJA,
961 Keeping in mind that the polynomial ``1`` evaluates the identity
962 element (also the zero element) of the trivial algebra, it is clear
963 that the polynomial ``1`` is the minimal polynomial of the only
964 element in a trivial algebra::
966 sage: J = TrivialEJA()
967 sage: J.one().minimal_polynomial()
969 sage: J.zero().minimal_polynomial()
974 The minimal polynomial of the identity and zero elements are
975 always the same, except in trivial algebras where the minimal
976 polynomial of the unit/zero element is ``1``::
978 sage: J = random_eja()
979 sage: mu = J.one().minimal_polynomial()
980 sage: t = mu.parent().gen()
981 sage: mu + int(J.is_trivial())*(t-2)
983 sage: mu = J.zero().minimal_polynomial()
984 sage: t = mu.parent().gen()
985 sage: mu + int(J.is_trivial())*(t-1)
988 The degree of an element is (by one definition) the degree
989 of its minimal polynomial::
991 sage: x = random_eja().random_element()
992 sage: x.degree() == x.minimal_polynomial().degree()
995 The minimal polynomial and the characteristic polynomial coincide
996 and are known (see Alizadeh, Example 11.11) for all elements of
997 the spin factor algebra that aren't scalar multiples of the
998 identity. We require the dimension of the algebra to be at least
999 two here so that said elements actually exist::
1001 sage: d_max = JordanSpinEJA._max_random_instance_dimension()
1002 sage: n = ZZ.random_element(2, max(2,d_max))
1003 sage: J = JordanSpinEJA(n)
1004 sage: y = J.random_element()
1005 sage: while y == y.coefficient(0)*J.one():
1006 ....: y = J.random_element()
1007 sage: y0 = y.to_vector()[0]
1008 sage: y_bar = y.to_vector()[1:]
1009 sage: actual = y.minimal_polynomial()
1010 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1011 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1012 sage: bool(actual == expected)
1015 The minimal polynomial should always kill its element::
1017 sage: x = random_eja().random_element() # long time
1018 sage: p = x.minimal_polynomial() # long time
1019 sage: x.apply_univariate_polynomial(p) # long time
1022 The minimal polynomial is invariant under a change of basis,
1023 and in particular, a re-scaling of the basis::
1025 sage: d_max = RealSymmetricEJA._max_random_instance_dimension()
1026 sage: d = ZZ.random_element(1, d_max)
1027 sage: n = RealSymmetricEJA._max_random_instance_size(d)
1028 sage: J1 = RealSymmetricEJA(n)
1029 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
1030 sage: X = random_matrix(AA,n)
1031 sage: X = X*X.transpose()
1034 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
1039 # Pretty sure we know what the minimal polynomial of
1040 # the zero operator is going to be. This ensures
1041 # consistency of e.g. the polynomial variable returned
1042 # in the "normal" case without us having to think about it.
1043 return self
.operator().minimal_polynomial()
1045 # If we don't orthonormalize the subalgebra's basis, then the
1046 # first two monomials in the subalgebra will be self^0 and
1047 # self^1... assuming that self^1 is not a scalar multiple of
1048 # self^0 (the unit element). We special case these to avoid
1049 # having to solve a system to coerce self into the subalgebra.
1050 A
= self
.subalgebra_generated_by(orthonormalize
=False)
1052 if A
.dimension() == 1:
1053 # Does a solve to find the scalar multiple alpha such that
1054 # alpha*unit = self. We have to do this because the basis
1055 # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
1056 unit
= self
.parent().one()
1057 alpha
= self
.to_vector() / unit
.to_vector()
1058 return (unit
.operator()*alpha
).minimal_polynomial()
1060 # If the dimension of the subalgebra is >= 2, then we just
1061 # use the second basis element.
1062 return A
.monomial(1).operator().minimal_polynomial()
1066 def to_matrix(self
):
1068 Return an (often more natural) representation of this element as a
1071 Every finite-dimensional Euclidean Jordan Algebra is a direct
1072 sum of five simple algebras, four of which comprise Hermitian
1073 matrices. This method returns a "natural" matrix
1074 representation of this element as either a Hermitian matrix or
1079 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1081 ....: QuaternionHermitianEJA,
1082 ....: RealSymmetricEJA)
1086 sage: J = ComplexHermitianEJA(3)
1089 sage: J.one().to_matrix()
1100 sage: J = QuaternionHermitianEJA(2)
1103 sage: J.one().to_matrix()
1110 This also works in Cartesian product algebras::
1112 sage: J1 = HadamardEJA(1)
1113 sage: J2 = RealSymmetricEJA(2)
1114 sage: J = cartesian_product([J1,J2])
1115 sage: x = sum(J.gens())
1116 sage: x.to_matrix()[0]
1118 sage: x.to_matrix()[1]
1119 [ 1 0.7071067811865475?]
1120 [0.7071067811865475? 1]
1123 B
= self
.parent().matrix_basis()
1124 W
= self
.parent().matrix_space()
1126 # This is just a manual "from_vector()", but of course
1127 # matrix spaces aren't vector spaces in sage, so they
1128 # don't have a from_vector() method.
1129 return W
.linear_combination( zip(B
, self
.to_vector()) )
1135 The norm of this element with respect to :meth:`inner_product`.
1139 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1144 sage: J = HadamardEJA(2)
1145 sage: x = sum(J.gens())
1151 sage: J = JordanSpinEJA(4)
1152 sage: x = sum(J.gens())
1157 return self
.inner_product(self
).sqrt()
1162 Return the left-multiplication-by-this-element
1163 operator on the ambient algebra.
1167 sage: from mjo.eja.eja_algebra import random_eja
1171 sage: J = random_eja()
1172 sage: x,y = J.random_elements(2)
1173 sage: x.operator()(y) == x*y
1175 sage: y.operator()(x) == x*y
1180 left_mult_by_self
= lambda y
: self
*y
1181 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1182 return EJAOperator(P
, P
, L
.matrix() )
1185 def quadratic_representation(self
, other
=None):
1187 Return the quadratic representation of this element.
1191 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1196 The explicit form in the spin factor algebra is given by
1197 Alizadeh's Example 11.12::
1199 sage: x = JordanSpinEJA.random_instance().random_element()
1200 sage: x_vec = x.to_vector()
1201 sage: Q = matrix.identity(x.base_ring(), 0)
1202 sage: n = x_vec.degree()
1205 ....: x_bar = x_vec[1:]
1206 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1207 ....: B = 2*x0*x_bar.row()
1208 ....: C = 2*x0*x_bar.column()
1209 ....: D = matrix.identity(x.base_ring(), n-1)
1210 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1211 ....: D = D + 2*x_bar.tensor_product(x_bar)
1212 ....: Q = matrix.block(2,2,[A,B,C,D])
1213 sage: Q == x.quadratic_representation().matrix()
1216 Test all of the properties from Theorem 11.2 in Alizadeh::
1218 sage: J = random_eja()
1219 sage: x,y = J.random_elements(2)
1220 sage: Lx = x.operator()
1221 sage: Lxx = (x*x).operator()
1222 sage: Qx = x.quadratic_representation()
1223 sage: Qy = y.quadratic_representation()
1224 sage: Qxy = x.quadratic_representation(y)
1225 sage: Qex = J.one().quadratic_representation(x)
1226 sage: n = ZZ.random_element(10)
1227 sage: Qxn = (x^n).quadratic_representation()
1231 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1234 Property 2 (multiply on the right for :trac:`28272`):
1236 sage: alpha = J.base_ring().random_element()
1237 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1242 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1245 sage: not x.is_invertible() or (
1248 ....: x.inverse().quadratic_representation() )
1251 sage: Qxy(J.one()) == x*y
1256 sage: not x.is_invertible() or (
1257 ....: x.quadratic_representation(x.inverse())*Qx
1258 ....: == Qx*x.quadratic_representation(x.inverse()) )
1261 sage: not x.is_invertible() or (
1262 ....: x.quadratic_representation(x.inverse())*Qx
1264 ....: 2*Lx*Qex - Qx )
1267 sage: 2*Lx*Qex - Qx == Lxx
1272 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1282 sage: not x.is_invertible() or (
1283 ....: Qx*x.inverse().operator() == Lx )
1288 sage: not x.operator_commutes_with(y) or (
1289 ....: Qx(y)^n == Qxn(y^n) )
1295 elif not other
in self
.parent():
1296 raise TypeError("'other' must live in the same algebra")
1299 M
= other
.operator()
1300 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1304 def spectral_decomposition(self
):
1306 Return the unique spectral decomposition of this element.
1310 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1311 element's left-multiplication-by operator to the subalgebra it
1312 generates. We then compute the spectral decomposition of that
1313 operator, and the spectral projectors we get back must be the
1314 left-multiplication-by operators for the idempotents we
1315 seek. Thus applying them to the identity element gives us those
1318 Since the eigenvalues are required to be distinct, we take
1319 the spectral decomposition of the zero element to be zero
1320 times the identity element of the algebra (which is idempotent,
1325 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1329 The spectral decomposition of the identity is ``1`` times itself,
1330 and the spectral decomposition of zero is ``0`` times the identity::
1332 sage: J = RealSymmetricEJA(3)
1335 sage: J.one().spectral_decomposition()
1337 sage: J.zero().spectral_decomposition()
1342 sage: J = RealSymmetricEJA(4)
1343 sage: x = sum(J.gens())
1344 sage: sd = x.spectral_decomposition()
1349 sage: c0.inner_product(c1) == 0
1351 sage: c0.is_idempotent()
1353 sage: c1.is_idempotent()
1355 sage: c0 + c1 == J.one()
1357 sage: l0*c0 + l1*c1 == x
1360 The spectral decomposition should work in subalgebras, too::
1362 sage: J = RealSymmetricEJA(4)
1363 sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
1364 sage: A = 2*b5 - 2*b8
1365 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1366 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1367 sage: (f0, f1, f2) = J1.gens()
1368 sage: f0.spectral_decomposition()
1369 [(0, 1.000000000000000?*c2), (1, 1.000000000000000?*c0)]
1372 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1374 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1375 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1378 def subalgebra_generated_by(self
, **kwargs
):
1380 Return the associative subalgebra of the parent EJA generated
1383 Since our parent algebra is unital, we want "subalgebra" to mean
1384 "unital subalgebra" as well; thus the subalgebra that an element
1385 generates will itself be a Euclidean Jordan algebra after
1386 restricting the algebra operations appropriately. This is the
1387 subalgebra that Faraut and Korányi work with in section II.2, for
1392 sage: from mjo.eja.eja_algebra import (random_eja,
1394 ....: RealSymmetricEJA)
1398 We can create subalgebras of Cartesian product EJAs that are not
1399 themselves Cartesian product EJAs (they're just "regular" EJAs)::
1401 sage: J1 = HadamardEJA(3)
1402 sage: J2 = RealSymmetricEJA(2)
1403 sage: J = cartesian_product([J1,J2])
1404 sage: J.one().subalgebra_generated_by()
1405 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
1409 This subalgebra, being composed of only powers, is associative::
1411 sage: x0 = random_eja().random_element()
1412 sage: A = x0.subalgebra_generated_by(orthonormalize=False)
1413 sage: x,y,z = A.random_elements(3)
1414 sage: (x*y)*z == x*(y*z)
1417 Squaring in the subalgebra should work the same as in
1420 sage: x = random_eja().random_element()
1421 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1422 sage: A(x^2) == A(x)*A(x)
1425 By definition, the subalgebra generated by the zero element is
1426 the one-dimensional algebra generated by the identity
1427 element... unless the original algebra was trivial, in which
1428 case the subalgebra is trivial too::
1430 sage: A = random_eja().zero().subalgebra_generated_by()
1431 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1435 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1436 A
= self
.parent().subalgebra(powers
,
1441 A
.one
.set_cache(A(self
.parent().one()))
1445 def subalgebra_idempotent(self
):
1447 Find an idempotent in the associative subalgebra I generate
1448 using Proposition 2.3.5 in Baes.
1452 sage: from mjo.eja.eja_algebra import random_eja
1456 Ensure that we can find an idempotent in a non-trivial algebra
1457 where there are non-nilpotent elements, or that we get the dumb
1458 solution in the trivial algebra::
1460 sage: J = random_eja(field=QQ, orthonormalize=False)
1461 sage: x = J.random_element()
1462 sage: while x.is_nilpotent() and not J.is_trivial():
1463 ....: x = J.random_element()
1464 sage: c = x.subalgebra_idempotent()
1469 if self
.parent().is_trivial():
1472 if self
.is_nilpotent():
1473 raise ValueError("this only works with non-nilpotent elements!")
1475 # The subalgebra is transient (we return an element of the
1476 # superalgebra, i.e. this algebra) so why bother
1478 J
= self
.subalgebra_generated_by(orthonormalize
=False)
1481 # The image of the matrix of left-u^m-multiplication
1482 # will be minimal for some natural number s...
1484 minimal_dim
= J
.dimension()
1485 for i
in range(1, minimal_dim
):
1486 this_dim
= (u
**i
).operator().matrix().image().dimension()
1487 if this_dim
< minimal_dim
:
1488 minimal_dim
= this_dim
1491 # Now minimal_matrix should correspond to the smallest
1492 # non-zero subspace in Baes's (or really, Koecher's)
1495 # However, we need to restrict the matrix to work on the
1496 # subspace... or do we? Can't we just solve, knowing that
1497 # A(c) = u^(s+1) should have a solution in the big space,
1500 A
= u_next
.operator().matrix()
1501 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1503 # Now c is the idempotent we want, but it still lives in
1505 return c
.superalgebra_element()
1510 Return my trace, the sum of my eigenvalues.
1512 In a trivial algebra, however you want to look at it, the trace is
1513 an empty sum for which we declare the result to be zero.
1517 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1524 sage: J = TrivialEJA()
1525 sage: J.zero().trace()
1529 sage: J = JordanSpinEJA(3)
1530 sage: x = sum(J.gens())
1536 sage: J = HadamardEJA(5)
1537 sage: J.one().trace()
1542 The trace of an element is a real number::
1544 sage: J = random_eja()
1545 sage: J.random_element().trace() in RLF
1548 The trace is linear::
1550 sage: J = random_eja()
1551 sage: x,y = J.random_elements(2)
1552 sage: alpha = J.base_ring().random_element()
1553 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1556 The trace of a square is nonnegative::
1558 sage: x = random_eja().random_element()
1559 sage: (x*x).trace() >= 0
1567 # Special case for the trivial algebra where
1568 # the trace is an empty sum.
1569 return P
.base_ring().zero()
1571 p
= P
._charpoly
_coefficients
()[r
-1]
1572 # The _charpoly_coeff function already adds the factor of
1573 # -1 to ensure that _charpoly_coeff(r-1) is really what
1574 # appears in front of t^{r-1} in the charpoly. However,
1575 # we want the negative of THAT for the trace.
1576 return -p(*self
.to_vector())
1579 def trace_inner_product(self
, other
):
1581 Return the trace inner product of myself and ``other``.
1585 sage: from mjo.eja.eja_algebra import random_eja
1589 The trace inner product is commutative, bilinear, and associative::
1591 sage: J = random_eja()
1592 sage: x,y,z = J.random_elements(3)
1594 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1597 sage: a = J.base_ring().random_element()
1598 sage: actual = (a*(x+z)).trace_inner_product(y)
1599 sage: expected = ( a*x.trace_inner_product(y) +
1600 ....: a*z.trace_inner_product(y) )
1601 sage: actual == expected
1603 sage: actual = x.trace_inner_product(a*(y+z))
1604 sage: expected = ( a*x.trace_inner_product(y) +
1605 ....: a*x.trace_inner_product(z) )
1606 sage: actual == expected
1609 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1613 if not other
in self
.parent():
1614 raise TypeError("'other' must live in the same algebra")
1616 return (self
*other
).trace()
1619 def trace_norm(self
):
1621 The norm of this element with respect to :meth:`trace_inner_product`.
1625 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1630 sage: J = HadamardEJA(2)
1631 sage: x = sum(J.gens())
1632 sage: x.trace_norm()
1637 sage: J = JordanSpinEJA(4)
1638 sage: x = sum(J.gens())
1639 sage: x.trace_norm()
1643 return self
.trace_inner_product(self
).sqrt()
1646 def operator_trace_inner_product(self
, other
):
1648 Return the operator inner product of myself and ``other``.
1650 The "operator inner product," whose name is not standard, is
1651 defined be the usual linear-algebraic trace of the
1652 ``(x*y).operator()``.
1654 Proposition III.1.5 in Faraut and Korányi shows that on any
1655 Euclidean Jordan algebra, this is another associative inner
1656 product under which the cone of squares is symmetric.
1658 This works even if the basis hasn't been orthonormalized
1659 because the eigenvalues of the corresponding matrix don't
1660 change when the basis does (they're preserved by any
1661 similarity transformation).
1665 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1666 ....: RealSymmetricEJA,
1667 ....: ComplexHermitianEJA,
1672 Proposition III.4.2 of Faraut and Korányi shows that on a
1673 simple algebra of rank `r` and dimension `n`, this inner
1674 product is `n/r` times the canonical
1675 :meth:`trace_inner_product`::
1677 sage: J = JordanSpinEJA(4, field=QQ)
1678 sage: x,y = J.random_elements(2)
1679 sage: n = J.dimension()
1681 sage: actual = x.operator_trace_inner_product(y)
1682 sage: expected = (n/r)*x.trace_inner_product(y)
1683 sage: actual == expected
1688 sage: J = RealSymmetricEJA(3)
1689 sage: x,y = J.random_elements(2)
1690 sage: n = J.dimension()
1692 sage: actual = x.operator_trace_inner_product(y)
1693 sage: expected = (n/r)*x.trace_inner_product(y)
1694 sage: actual == expected
1699 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
1700 sage: x,y = J.random_elements(2)
1701 sage: n = J.dimension()
1703 sage: actual = x.operator_trace_inner_product(y)
1704 sage: expected = (n/r)*x.trace_inner_product(y)
1705 sage: actual == expected
1710 The operator inner product is commutative, bilinear, and
1713 sage: J = random_eja()
1714 sage: x,y,z = J.random_elements(3)
1716 sage: actual = x.operator_trace_inner_product(y)
1717 sage: expected = y.operator_trace_inner_product(x)
1718 sage: actual == expected
1721 sage: a = J.base_ring().random_element()
1722 sage: actual = (a*(x+z)).operator_trace_inner_product(y)
1723 sage: expected = ( a*x.operator_trace_inner_product(y) +
1724 ....: a*z.operator_trace_inner_product(y) )
1725 sage: actual == expected
1727 sage: actual = x.operator_trace_inner_product(a*(y+z))
1728 sage: expected = ( a*x.operator_trace_inner_product(y) +
1729 ....: a*x.operator_trace_inner_product(z) )
1730 sage: actual == expected
1733 sage: actual = (x*y).operator_trace_inner_product(z)
1734 sage: expected = y.operator_trace_inner_product(x*z)
1735 sage: actual == expected
1738 Despite the fact that the implementation uses a matrix representation,
1739 the answer is independent of the basis used::
1741 sage: J = RealSymmetricEJA(3, field=QQ, orthonormalize=False)
1742 sage: V = RealSymmetricEJA(3)
1743 sage: x,y = J.random_elements(2)
1744 sage: w = V(x.to_matrix())
1745 sage: z = V(y.to_matrix())
1746 sage: expected = x.operator_trace_inner_product(y)
1747 sage: actual = w.operator_trace_inner_product(z)
1748 sage: actual == expected
1752 if not other
in self
.parent():
1753 raise TypeError("'other' must live in the same algebra")
1755 return (self
*other
).operator().matrix().trace()
1758 def operator_trace_norm(self
):
1760 The norm of this element with respect to
1761 :meth:`operator_trace_inner_product`.
1765 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1770 On a simple algebra, this will differ from :meth:`trace_norm`
1771 by the scalar factor ``(n/r).sqrt()``, where `n` is the
1772 dimension of the algebra and `r` its rank. This follows from
1773 the corresponding result (Proposition III.4.2 of Faraut and
1774 Korányi) for the trace inner product::
1776 sage: J = HadamardEJA(2)
1777 sage: x = sum(J.gens())
1778 sage: x.operator_trace_norm()
1783 sage: J = JordanSpinEJA(4)
1784 sage: x = sum(J.gens())
1785 sage: x.operator_trace_norm()
1789 return self
.operator_trace_inner_product(self
).sqrt()
1792 class CartesianProductParentEJAElement(EJAElement
):
1794 An intermediate class for elements that have a Cartesian
1795 product as their parent algebra.
1797 This is needed because the ``to_matrix`` method (which gives you a
1798 representation from the superalgebra) needs to do special stuff
1799 for Cartesian products. Specifically, an EJA subalgebra of a
1800 Cartesian product EJA will not itself be a Cartesian product (it
1801 has its own basis) -- but we want ``to_matrix()`` to be able to
1802 give us a Cartesian product representation.
1804 def to_matrix(self
):
1805 # An override is necessary to call our custom _scale().
1806 B
= self
.parent().matrix_basis()
1807 W
= self
.parent().matrix_space()
1809 # Aaaaand linear combinations don't work in Cartesian
1810 # product spaces, even though they provide a method with
1811 # that name. This is hidden in a subclass because the
1812 # _scale() function is slow.
1813 pairs
= zip(B
, self
.to_vector())
1814 return W
.sum( _scale(b
, alpha
) for (b
,alpha
) in pairs
)
1816 class CartesianProductEJAElement(CartesianProductParentEJAElement
):
1819 Compute the determinant of this product-element using the
1820 determianants of its factors.
1822 This result Follows from the spectral decomposition of (say)
1823 the pair `(x,y)` in terms of the Jordan frame `\left\{ (c_1,
1824 0),(c_2, 0),...,(0,d_1),(0,d_2),... \right\}.
1826 from sage
.misc
.misc_c
import prod
1827 return prod( f
.det() for f
in self
.cartesian_factors() )