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 _scale
10 class FiniteDimensionalEJAElement(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: set_random_seed()
47 sage: x = random_eja().random_element()
51 A few examples of power-associativity::
53 sage: set_random_seed()
54 sage: x = random_eja().random_element()
55 sage: x*(x*x)*(x*x) == x^5
57 sage: (x*x)*(x*x*x) == x^5
60 We also know that powers operator-commute (Koecher, Chapter
63 sage: set_random_seed()
64 sage: x = random_eja().random_element()
65 sage: m = ZZ.random_element(0,10)
66 sage: n = ZZ.random_element(0,10)
67 sage: Lxm = (x^m).operator()
68 sage: Lxn = (x^n).operator()
69 sage: Lxm*Lxn == Lxn*Lxm
74 return self
.parent().one()
78 return (self
**(n
-1))*self
81 def apply_univariate_polynomial(self
, p
):
83 Apply the univariate polynomial ``p`` to this element.
85 A priori, SageMath won't allow us to apply a univariate
86 polynomial to an element of an EJA, because we don't know
87 that EJAs are rings (they are usually not associative). Of
88 course, we know that EJAs are power-associative, so the
89 operation is ultimately kosher. This function sidesteps
90 the CAS to get the answer we want and expect.
94 sage: from mjo.eja.eja_algebra import (HadamardEJA,
99 sage: R = PolynomialRing(QQ, 't')
101 sage: p = t^4 - t^3 + 5*t - 2
102 sage: J = HadamardEJA(5)
103 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
108 We should always get back an element of the algebra::
110 sage: set_random_seed()
111 sage: p = PolynomialRing(AA, 't').random_element()
112 sage: J = random_eja()
113 sage: x = J.random_element()
114 sage: x.apply_univariate_polynomial(p) in J
118 if len(p
.variables()) > 1:
119 raise ValueError("not a univariate polynomial")
122 # Convert the coeficcients to the parent's base ring,
123 # because a priori they might live in an (unnecessarily)
124 # larger ring for which P.sum() would fail below.
125 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
126 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
129 def characteristic_polynomial(self
):
131 Return the characteristic polynomial of this element.
135 sage: from mjo.eja.eja_algebra import (random_eja,
140 The rank of `R^3` is three, and the minimal polynomial of
141 the identity element is `(t-1)` from which it follows that
142 the characteristic polynomial should be `(t-1)^3`::
144 sage: J = HadamardEJA(3)
145 sage: J.one().characteristic_polynomial()
146 t^3 - 3*t^2 + 3*t - 1
148 Likewise, the characteristic of the zero element in the
149 rank-three algebra `R^{n}` should be `t^{3}`::
151 sage: J = HadamardEJA(3)
152 sage: J.zero().characteristic_polynomial()
157 The characteristic polynomial of an element should evaluate
158 to zero on that element::
160 sage: set_random_seed()
161 sage: x = random_eja().random_element()
162 sage: p = x.characteristic_polynomial()
163 sage: x.apply_univariate_polynomial(p).is_zero()
166 The characteristic polynomials of the zero and unit elements
167 should be what we think they are in a subalgebra, too::
169 sage: J = HadamardEJA(3)
170 sage: p1 = J.one().characteristic_polynomial()
171 sage: q1 = J.zero().characteristic_polynomial()
172 sage: b0,b1,b2 = J.gens()
173 sage: A = (b0 + 2*b1 + 3*b2).subalgebra_generated_by() # dim 3
174 sage: p2 = A.one().characteristic_polynomial()
175 sage: q2 = A.zero().characteristic_polynomial()
182 p
= self
.parent().characteristic_polynomial_of()
183 return p(*self
.to_vector())
186 def inner_product(self
, other
):
188 Return the parent algebra's inner product of myself and ``other``.
192 sage: from mjo.eja.eja_algebra import (
193 ....: ComplexHermitianEJA,
195 ....: QuaternionHermitianEJA,
196 ....: RealSymmetricEJA,
201 The inner product in the Jordan spin algebra is the usual
202 inner product on `R^n` (this example only works because the
203 basis for the Jordan algebra is the standard basis in `R^n`)::
205 sage: J = JordanSpinEJA(3)
206 sage: x = vector(QQ,[1,2,3])
207 sage: y = vector(QQ,[4,5,6])
208 sage: x.inner_product(y)
210 sage: J.from_vector(x).inner_product(J.from_vector(y))
213 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
214 multiplication is the usual matrix multiplication in `S^n`,
215 so the inner product of the identity matrix with itself
218 sage: J = RealSymmetricEJA(3)
219 sage: J.one().inner_product(J.one())
222 Likewise, the inner product on `C^n` is `<X,Y> =
223 Re(trace(X*Y))`, where we must necessarily take the real
224 part because the product of Hermitian matrices may not be
227 sage: J = ComplexHermitianEJA(3)
228 sage: J.one().inner_product(J.one())
231 Ditto for the quaternions::
233 sage: J = QuaternionHermitianEJA(2)
234 sage: J.one().inner_product(J.one())
239 Ensure that we can always compute an inner product, and that
240 it gives us back a real number::
242 sage: set_random_seed()
243 sage: J = random_eja()
244 sage: x,y = J.random_elements(2)
245 sage: x.inner_product(y) in RLF
251 raise TypeError("'other' must live in the same algebra")
253 return P
.inner_product(self
, other
)
256 def operator_commutes_with(self
, other
):
258 Return whether or not this element operator-commutes
263 sage: from mjo.eja.eja_algebra import random_eja
267 The definition of a Jordan algebra says that any element
268 operator-commutes with its square::
270 sage: set_random_seed()
271 sage: x = random_eja().random_element()
272 sage: x.operator_commutes_with(x^2)
277 Test Lemma 1 from Chapter III of Koecher::
279 sage: set_random_seed()
280 sage: u,v = random_eja().random_elements(2)
281 sage: lhs = u.operator_commutes_with(u*v)
282 sage: rhs = v.operator_commutes_with(u^2)
286 Test the first polarization identity from my notes, Koecher
287 Chapter III, or from Baes (2.3)::
289 sage: set_random_seed()
290 sage: x,y = random_eja().random_elements(2)
291 sage: Lx = x.operator()
292 sage: Ly = y.operator()
293 sage: Lxx = (x*x).operator()
294 sage: Lxy = (x*y).operator()
295 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
298 Test the second polarization identity from my notes or from
301 sage: set_random_seed()
302 sage: x,y,z = random_eja().random_elements(3)
303 sage: Lx = x.operator()
304 sage: Ly = y.operator()
305 sage: Lz = z.operator()
306 sage: Lzy = (z*y).operator()
307 sage: Lxy = (x*y).operator()
308 sage: Lxz = (x*z).operator()
309 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
312 Test the third polarization identity from my notes or from
315 sage: set_random_seed()
316 sage: u,y,z = random_eja().random_elements(3)
317 sage: Lu = u.operator()
318 sage: Ly = y.operator()
319 sage: Lz = z.operator()
320 sage: Lzy = (z*y).operator()
321 sage: Luy = (u*y).operator()
322 sage: Luz = (u*z).operator()
323 sage: Luyz = (u*(y*z)).operator()
324 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
325 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
326 sage: bool(lhs == rhs)
330 if not other
in self
.parent():
331 raise TypeError("'other' must live in the same algebra")
340 Return my determinant, the product of my eigenvalues.
344 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
346 ....: RealSymmetricEJA,
347 ....: ComplexHermitianEJA,
352 sage: J = JordanSpinEJA(2)
353 sage: x = sum( J.gens() )
359 sage: J = JordanSpinEJA(3)
360 sage: x = sum( J.gens() )
364 The determinant of the sole element in the rank-zero trivial
365 algebra is ``1``, by three paths of reasoning. First, its
366 characteristic polynomial is a constant ``1``, so the constant
367 term in that polynomial is ``1``. Second, the characteristic
368 polynomial evaluated at zero is again ``1``. And finally, the
369 (empty) product of its eigenvalues is likewise just unity::
371 sage: J = TrivialEJA()
377 An element is invertible if and only if its determinant is
380 sage: set_random_seed()
381 sage: x = random_eja().random_element()
382 sage: x.is_invertible() == (x.det() != 0)
385 Ensure that the determinant is multiplicative on an associative
386 subalgebra as in Faraut and Korányi's Proposition II.2.2::
388 sage: set_random_seed()
389 sage: J = random_eja().random_element().subalgebra_generated_by()
390 sage: x,y = J.random_elements(2)
391 sage: (x*y).det() == x.det()*y.det()
394 The determinant in real matrix algebras is the usual determinant::
396 sage: set_random_seed()
397 sage: X = matrix.random(QQ,3)
399 sage: J1 = RealSymmetricEJA(3)
400 sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
401 sage: expected = X.det()
402 sage: actual1 = J1(X).det()
403 sage: actual2 = J2(X).det()
404 sage: actual1 == expected
406 sage: actual2 == expected
414 # Special case, since we don't get the a0=1
415 # coefficient when the rank of the algebra
417 return P
.base_ring().one()
419 p
= P
._charpoly
_coefficients
()[0]
420 # The _charpoly_coeff function already adds the factor of -1
421 # to ensure that _charpoly_coefficients()[0] is really what
422 # appears in front of t^{0} in the charpoly. However, we want
423 # (-1)^r times THAT for the determinant.
424 return ((-1)**r
)*p(*self
.to_vector())
430 Return the Jordan-multiplicative inverse of this element.
434 In general we appeal to the quadratic representation as in
435 Koecher's Theorem 12 in Chapter III, Section 5. But if the
436 parent algebra's "characteristic polynomial of" coefficients
437 happen to be cached, then we use Proposition II.2.4 in Faraut
438 and Korányi which gives a formula for the inverse based on the
439 characteristic polynomial and the Cayley-Hamilton theorem for
440 Euclidean Jordan algebras::
444 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
450 The inverse in the spin factor algebra is given in Alizadeh's
453 sage: set_random_seed()
454 sage: J = JordanSpinEJA.random_instance()
455 sage: x = J.random_element()
456 sage: while not x.is_invertible():
457 ....: x = J.random_element()
458 sage: x_vec = x.to_vector()
460 sage: x_bar = x_vec[1:]
461 sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
462 sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
463 sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
464 sage: x.inverse() == J.from_vector(x_inverse)
467 Trying to invert a non-invertible element throws an error:
469 sage: JordanSpinEJA(3).zero().inverse()
470 Traceback (most recent call last):
472 ZeroDivisionError: element is not invertible
476 The identity element is its own inverse::
478 sage: set_random_seed()
479 sage: J = random_eja()
480 sage: J.one().inverse() == J.one()
483 If an element has an inverse, it acts like one::
485 sage: set_random_seed()
486 sage: J = random_eja()
487 sage: x = J.random_element()
488 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
491 The inverse of the inverse is what we started with::
493 sage: set_random_seed()
494 sage: J = random_eja()
495 sage: x = J.random_element()
496 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
499 Proposition II.2.3 in Faraut and Korányi says that the inverse
500 of an element is the inverse of its left-multiplication operator
501 applied to the algebra's identity, when that inverse exists::
503 sage: set_random_seed()
504 sage: J = random_eja()
505 sage: x = J.random_element()
506 sage: (not x.operator().is_invertible()) or (
507 ....: x.operator().inverse()(J.one()) == x.inverse() )
510 Check that the fast (cached) and slow algorithms give the same
513 sage: set_random_seed() # long time
514 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
515 sage: x = J.random_element() # long time
516 sage: while not x.is_invertible(): # long time
517 ....: x = J.random_element() # long time
518 sage: slow = x.inverse() # long time
519 sage: _ = J._charpoly_coefficients() # long time
520 sage: fast = x.inverse() # long time
521 sage: slow == fast # long time
524 not_invertible_msg
= "element is not invertible"
525 if self
.parent()._charpoly
_coefficients
.is_in_cache():
526 # We can invert using our charpoly if it will be fast to
527 # compute. If the coefficients are cached, our rank had
529 if self
.det().is_zero():
530 raise ZeroDivisionError(not_invertible_msg
)
531 r
= self
.parent().rank()
532 a
= self
.characteristic_polynomial().coefficients(sparse
=False)
533 return (-1)**(r
+1)*sum(a
[i
+1]*self
**i
for i
in range(r
))/self
.det()
536 inv
= (~self
.quadratic_representation())(self
)
537 self
.is_invertible
.set_cache(True)
539 except ZeroDivisionError:
540 self
.is_invertible
.set_cache(False)
541 raise ZeroDivisionError(not_invertible_msg
)
545 def is_invertible(self
):
547 Return whether or not this element is invertible.
551 If computing my determinant will be fast, we do so and compare
552 with zero (Proposition II.2.4 in Faraut and
553 Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
554 reduces the problem to the invertibility of my quadratic
559 sage: from mjo.eja.eja_algebra import random_eja
563 The identity element is always invertible::
565 sage: set_random_seed()
566 sage: J = random_eja()
567 sage: J.one().is_invertible()
570 The zero element is never invertible in a non-trivial algebra::
572 sage: set_random_seed()
573 sage: J = random_eja()
574 sage: (not J.is_trivial()) and J.zero().is_invertible()
577 Test that the fast (cached) and slow algorithms give the same
580 sage: set_random_seed() # long time
581 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
582 sage: x = J.random_element() # long time
583 sage: slow = x.is_invertible() # long time
584 sage: _ = J._charpoly_coefficients() # long time
585 sage: fast = x.is_invertible() # long time
586 sage: slow == fast # long time
590 if self
.parent().is_trivial():
595 if self
.parent()._charpoly
_coefficients
.is_in_cache():
596 # The determinant will be quicker than inverting the
597 # quadratic representation, most likely.
598 return (not self
.det().is_zero())
600 # The easiest way to determine if I'm invertible is to try.
602 inv
= (~self
.quadratic_representation())(self
)
603 self
.inverse
.set_cache(inv
)
605 except ZeroDivisionError:
609 def is_primitive_idempotent(self
):
611 Return whether or not this element is a primitive (or minimal)
614 A primitive idempotent is a non-zero idempotent that is not
615 the sum of two other non-zero idempotents. Remark 2.7.15 in
616 Baes shows that this is what he refers to as a "minimal
619 An element of a Euclidean Jordan algebra is a minimal idempotent
620 if it :meth:`is_idempotent` and if its Peirce subalgebra
621 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
626 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
627 ....: RealSymmetricEJA,
633 This method is sloooooow.
637 The spectral decomposition of a non-regular element should always
638 contain at least one non-minimal idempotent::
640 sage: J = RealSymmetricEJA(3)
641 sage: x = sum(J.gens())
644 sage: [ c.is_primitive_idempotent()
645 ....: for (l,c) in x.spectral_decomposition() ]
648 On the other hand, the spectral decomposition of a regular
649 element should always be in terms of minimal idempotents::
651 sage: J = JordanSpinEJA(4)
652 sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
655 sage: [ c.is_primitive_idempotent()
656 ....: for (l,c) in x.spectral_decomposition() ]
661 The identity element is minimal only in an EJA of rank one::
663 sage: set_random_seed()
664 sage: J = random_eja()
665 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
668 A non-idempotent cannot be a minimal idempotent::
670 sage: set_random_seed()
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: set_random_seed()
681 sage: J = JordanSpinEJA(4)
682 sage: x = J.random_element()
683 sage: expected = (x.is_idempotent() and x.trace() == 1)
684 sage: actual = x.is_primitive_idempotent()
685 sage: actual == expected
688 Primitive idempotents must be non-zero::
690 sage: set_random_seed()
691 sage: J = random_eja()
692 sage: J.zero().is_idempotent()
694 sage: J.zero().is_primitive_idempotent()
697 As a consequence of the fact that primitive idempotents must
698 be non-zero, there are no primitive idempotents in a trivial
699 Euclidean Jordan algebra::
701 sage: J = TrivialEJA()
702 sage: J.one().is_idempotent()
704 sage: J.one().is_primitive_idempotent()
708 if not self
.is_idempotent():
714 (_
,_
,J1
) = self
.parent().peirce_decomposition(self
)
715 return (J1
.dimension() == 1)
718 def is_nilpotent(self
):
720 Return whether or not some power of this element is zero.
724 We use Theorem 5 in Chapter III of Koecher, which says that
725 an element ``x`` is nilpotent if and only if ``x.operator()``
726 is nilpotent. And it is a basic fact of linear algebra that
727 an operator on an `n`-dimensional space is nilpotent if and
728 only if, when raised to the `n`th power, it equals the zero
729 operator (for example, see Axler Corollary 8.8).
733 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
738 sage: J = JordanSpinEJA(3)
739 sage: x = sum(J.gens())
740 sage: x.is_nilpotent()
745 The identity element is never nilpotent, except in a trivial EJA::
747 sage: set_random_seed()
748 sage: J = random_eja()
749 sage: J.one().is_nilpotent() and not J.is_trivial()
752 The additive identity is always nilpotent::
754 sage: set_random_seed()
755 sage: random_eja().zero().is_nilpotent()
760 zero_operator
= P
.zero().operator()
761 return self
.operator()**P
.dimension() == zero_operator
764 def is_regular(self
):
766 Return whether or not this is a regular element.
770 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
775 The identity element always has degree one, but any element
776 linearly-independent from it is regular::
778 sage: J = JordanSpinEJA(5)
779 sage: J.one().is_regular()
781 sage: b0, b1, b2, b3, b4 = J.gens()
784 sage: for x in J.gens():
785 ....: (J.one() + x).is_regular()
794 The zero element should never be regular, unless the parent
795 algebra has dimension less than or equal to one::
797 sage: set_random_seed()
798 sage: J = random_eja()
799 sage: J.dimension() <= 1 or not J.zero().is_regular()
802 The unit element isn't regular unless the algebra happens to
803 consist of only its scalar multiples::
805 sage: set_random_seed()
806 sage: J = random_eja()
807 sage: J.dimension() <= 1 or not J.one().is_regular()
811 return self
.degree() == self
.parent().rank()
816 Return the degree of this element, which is defined to be
817 the degree of its minimal polynomial.
825 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
830 sage: J = JordanSpinEJA(4)
831 sage: J.one().degree()
833 sage: b0,b1,b2,b3 = J.gens()
834 sage: (b0 - b1).degree()
837 In the spin factor algebra (of rank two), all elements that
838 aren't multiples of the identity are regular::
840 sage: set_random_seed()
841 sage: J = JordanSpinEJA.random_instance()
842 sage: n = J.dimension()
843 sage: x = J.random_element()
844 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
849 The zero and unit elements are both of degree one in nontrivial
852 sage: set_random_seed()
853 sage: J = random_eja()
854 sage: d = J.zero().degree()
855 sage: (J.is_trivial() and d == 0) or d == 1
857 sage: d = J.one().degree()
858 sage: (J.is_trivial() and d == 0) or d == 1
861 Our implementation agrees with the definition::
863 sage: set_random_seed()
864 sage: x = random_eja().random_element()
865 sage: x.degree() == x.minimal_polynomial().degree()
869 n
= self
.parent().dimension()
872 # The minimal polynomial is an empty product, i.e. the
873 # constant polynomial "1" having degree zero.
876 # The minimal polynomial of zero in a nontrivial algebra
877 # is "t", and is of degree one.
880 # If this is a nonzero element of a nontrivial algebra, it
881 # has degree at least one. It follows that, in an algebra
882 # of dimension one, the degree must be actually one.
885 # BEWARE: The subalgebra_generated_by() method uses the result
886 # of this method to construct a basis for the subalgebra. That
887 # means, in particular, that we cannot implement this method
888 # as ``self.subalgebra_generated_by().dimension()``.
890 # Algorithm: keep appending (vector representations of) powers
891 # self as rows to a matrix and echelonizing it. When its rank
892 # stops increasing, we've reached a redundancy.
894 # Given the special cases above, we can assume that "self" is
895 # nonzero, the algebra is nontrivial, and that its dimension
897 M
= matrix([(self
.parent().one()).to_vector()])
900 # Specifying the row-reduction algorithm can e.g. help over
901 # AA because it avoids the RecursionError that gets thrown
902 # when we have to look too hard for a root.
904 # Beware: QQ supports an entirely different set of "algorithm"
905 # keywords than do AA and RR.
907 from sage
.rings
.all
import QQ
908 if self
.parent().base_ring() is not QQ
:
909 algo
= "scaled_partial_pivoting"
912 M
= matrix(M
.rows() + [(self
**d
).to_vector()])
915 if new_rank
== old_rank
:
924 def left_matrix(self
):
926 Our parent class defines ``left_matrix`` and ``matrix``
927 methods whose names are misleading. We don't want them.
929 raise NotImplementedError("use operator().matrix() instead")
934 def minimal_polynomial(self
):
936 Return the minimal polynomial of this element,
937 as a function of the variable `t`.
941 We restrict ourselves to the associative subalgebra
942 generated by this element, and then return the minimal
943 polynomial of this element's operator matrix (in that
944 subalgebra). This works by Baes Proposition 2.3.16.
948 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
949 ....: RealSymmetricEJA,
955 Keeping in mind that the polynomial ``1`` evaluates the identity
956 element (also the zero element) of the trivial algebra, it is clear
957 that the polynomial ``1`` is the minimal polynomial of the only
958 element in a trivial algebra::
960 sage: J = TrivialEJA()
961 sage: J.one().minimal_polynomial()
963 sage: J.zero().minimal_polynomial()
968 The minimal polynomial of the identity and zero elements are
969 always the same, except in trivial algebras where the minimal
970 polynomial of the unit/zero element is ``1``::
972 sage: set_random_seed()
973 sage: J = random_eja()
974 sage: mu = J.one().minimal_polynomial()
975 sage: t = mu.parent().gen()
976 sage: mu + int(J.is_trivial())*(t-2)
978 sage: mu = J.zero().minimal_polynomial()
979 sage: t = mu.parent().gen()
980 sage: mu + int(J.is_trivial())*(t-1)
983 The degree of an element is (by one definition) the degree
984 of its minimal polynomial::
986 sage: set_random_seed()
987 sage: x = random_eja().random_element()
988 sage: x.degree() == x.minimal_polynomial().degree()
991 The minimal polynomial and the characteristic polynomial coincide
992 and are known (see Alizadeh, Example 11.11) for all elements of
993 the spin factor algebra that aren't scalar multiples of the
994 identity. We require the dimension of the algebra to be at least
995 two here so that said elements actually exist::
997 sage: set_random_seed()
998 sage: d_max = JordanSpinEJA._max_random_instance_dimension()
999 sage: n = ZZ.random_element(2, max(2,d_max))
1000 sage: J = JordanSpinEJA(n)
1001 sage: y = J.random_element()
1002 sage: while y == y.coefficient(0)*J.one():
1003 ....: y = J.random_element()
1004 sage: y0 = y.to_vector()[0]
1005 sage: y_bar = y.to_vector()[1:]
1006 sage: actual = y.minimal_polynomial()
1007 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1008 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1009 sage: bool(actual == expected)
1012 The minimal polynomial should always kill its element::
1014 sage: set_random_seed()
1015 sage: x = random_eja().random_element()
1016 sage: p = x.minimal_polynomial()
1017 sage: x.apply_univariate_polynomial(p)
1020 The minimal polynomial is invariant under a change of basis,
1021 and in particular, a re-scaling of the basis::
1023 sage: set_random_seed()
1024 sage: d_max = RealSymmetricEJA._max_random_instance_dimension()
1025 sage: d = ZZ.random_element(1, d_max)
1026 sage: n = RealSymmetricEJA._max_random_instance_size(d)
1027 sage: J1 = RealSymmetricEJA(n)
1028 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
1029 sage: X = random_matrix(AA,n)
1030 sage: X = X*X.transpose()
1033 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
1038 # Pretty sure we know what the minimal polynomial of
1039 # the zero operator is going to be. This ensures
1040 # consistency of e.g. the polynomial variable returned
1041 # in the "normal" case without us having to think about it.
1042 return self
.operator().minimal_polynomial()
1044 # If we don't orthonormalize the subalgebra's basis, then the
1045 # first two monomials in the subalgebra will be self^0 and
1046 # self^1... assuming that self^1 is not a scalar multiple of
1047 # self^0 (the unit element). We special case these to avoid
1048 # having to solve a system to coerce self into the subalgebra.
1049 A
= self
.subalgebra_generated_by(orthonormalize
=False)
1051 if A
.dimension() == 1:
1052 # Does a solve to find the scalar multiple alpha such that
1053 # alpha*unit = self. We have to do this because the basis
1054 # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
1055 unit
= self
.parent().one()
1056 alpha
= self
.to_vector() / unit
.to_vector()
1057 return (unit
.operator()*alpha
).minimal_polynomial()
1059 # If the dimension of the subalgebra is >= 2, then we just
1060 # use the second basis element.
1061 return A
.monomial(1).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,
1080 ....: QuaternionHermitianEJA,
1081 ....: RealSymmetricEJA)
1085 sage: J = ComplexHermitianEJA(3)
1088 sage: J.one().to_matrix()
1099 sage: J = QuaternionHermitianEJA(2)
1102 sage: J.one().to_matrix()
1109 This also works in Cartesian product algebras::
1111 sage: J1 = HadamardEJA(1)
1112 sage: J2 = RealSymmetricEJA(2)
1113 sage: J = cartesian_product([J1,J2])
1114 sage: x = sum(J.gens())
1115 sage: x.to_matrix()[0]
1117 sage: x.to_matrix()[1]
1118 [ 1 0.7071067811865475?]
1119 [0.7071067811865475? 1]
1122 B
= self
.parent().matrix_basis()
1123 W
= self
.parent().matrix_space()
1125 # This is just a manual "from_vector()", but of course
1126 # matrix spaces aren't vector spaces in sage, so they
1127 # don't have a from_vector() method.
1128 return W
.linear_combination( zip(B
, self
.to_vector()) )
1134 The norm of this element with respect to :meth:`inner_product`.
1138 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1143 sage: J = HadamardEJA(2)
1144 sage: x = sum(J.gens())
1150 sage: J = JordanSpinEJA(4)
1151 sage: x = sum(J.gens())
1156 return self
.inner_product(self
).sqrt()
1161 Return the left-multiplication-by-this-element
1162 operator on the ambient algebra.
1166 sage: from mjo.eja.eja_algebra import random_eja
1170 sage: set_random_seed()
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 FiniteDimensionalEJAOperator(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: set_random_seed()
1200 sage: x = JordanSpinEJA.random_instance().random_element()
1201 sage: x_vec = x.to_vector()
1202 sage: Q = matrix.identity(x.base_ring(), 0)
1203 sage: n = x_vec.degree()
1206 ....: x_bar = x_vec[1:]
1207 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1208 ....: B = 2*x0*x_bar.row()
1209 ....: C = 2*x0*x_bar.column()
1210 ....: D = matrix.identity(x.base_ring(), n-1)
1211 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1212 ....: D = D + 2*x_bar.tensor_product(x_bar)
1213 ....: Q = matrix.block(2,2,[A,B,C,D])
1214 sage: Q == x.quadratic_representation().matrix()
1217 Test all of the properties from Theorem 11.2 in Alizadeh::
1219 sage: set_random_seed()
1220 sage: J = random_eja()
1221 sage: x,y = J.random_elements(2)
1222 sage: Lx = x.operator()
1223 sage: Lxx = (x*x).operator()
1224 sage: Qx = x.quadratic_representation()
1225 sage: Qy = y.quadratic_representation()
1226 sage: Qxy = x.quadratic_representation(y)
1227 sage: Qex = J.one().quadratic_representation(x)
1228 sage: n = ZZ.random_element(10)
1229 sage: Qxn = (x^n).quadratic_representation()
1233 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1236 Property 2 (multiply on the right for :trac:`28272`):
1238 sage: alpha = J.base_ring().random_element()
1239 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1244 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1247 sage: not x.is_invertible() or (
1250 ....: x.inverse().quadratic_representation() )
1253 sage: Qxy(J.one()) == x*y
1258 sage: not x.is_invertible() or (
1259 ....: x.quadratic_representation(x.inverse())*Qx
1260 ....: == Qx*x.quadratic_representation(x.inverse()) )
1263 sage: not x.is_invertible() or (
1264 ....: x.quadratic_representation(x.inverse())*Qx
1266 ....: 2*Lx*Qex - Qx )
1269 sage: 2*Lx*Qex - Qx == Lxx
1274 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1284 sage: not x.is_invertible() or (
1285 ....: Qx*x.inverse().operator() == Lx )
1290 sage: not x.operator_commutes_with(y) or (
1291 ....: Qx(y)^n == Qxn(y^n) )
1297 elif not other
in self
.parent():
1298 raise TypeError("'other' must live in the same algebra")
1301 M
= other
.operator()
1302 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1306 def spectral_decomposition(self
):
1308 Return the unique spectral decomposition of this element.
1312 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1313 element's left-multiplication-by operator to the subalgebra it
1314 generates. We then compute the spectral decomposition of that
1315 operator, and the spectral projectors we get back must be the
1316 left-multiplication-by operators for the idempotents we
1317 seek. Thus applying them to the identity element gives us those
1320 Since the eigenvalues are required to be distinct, we take
1321 the spectral decomposition of the zero element to be zero
1322 times the identity element of the algebra (which is idempotent,
1327 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1331 The spectral decomposition of the identity is ``1`` times itself,
1332 and the spectral decomposition of zero is ``0`` times the identity::
1334 sage: J = RealSymmetricEJA(3)
1337 sage: J.one().spectral_decomposition()
1339 sage: J.zero().spectral_decomposition()
1344 sage: J = RealSymmetricEJA(4)
1345 sage: x = sum(J.gens())
1346 sage: sd = x.spectral_decomposition()
1351 sage: c0.inner_product(c1) == 0
1353 sage: c0.is_idempotent()
1355 sage: c1.is_idempotent()
1357 sage: c0 + c1 == J.one()
1359 sage: l0*c0 + l1*c1 == x
1362 The spectral decomposition should work in subalgebras, too::
1364 sage: J = RealSymmetricEJA(4)
1365 sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
1366 sage: A = 2*b5 - 2*b8
1367 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1368 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1369 sage: (f0, f1, f2) = J1.gens()
1370 sage: f0.spectral_decomposition()
1371 [(0, 1.000000000000000?*c2), (1, 1.000000000000000?*c0)]
1374 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1376 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1377 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1380 def subalgebra_generated_by(self
, **kwargs
):
1382 Return the associative subalgebra of the parent EJA generated
1385 Since our parent algebra is unital, we want "subalgebra" to mean
1386 "unital subalgebra" as well; thus the subalgebra that an element
1387 generates will itself be a Euclidean Jordan algebra after
1388 restricting the algebra operations appropriately. This is the
1389 subalgebra that Faraut and Korányi work with in section II.2, for
1394 sage: from mjo.eja.eja_algebra import (random_eja,
1396 ....: RealSymmetricEJA)
1400 We can create subalgebras of Cartesian product EJAs that are not
1401 themselves Cartesian product EJAs (they're just "regular" EJAs)::
1403 sage: J1 = HadamardEJA(3)
1404 sage: J2 = RealSymmetricEJA(2)
1405 sage: J = cartesian_product([J1,J2])
1406 sage: J.one().subalgebra_generated_by()
1407 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
1411 This subalgebra, being composed of only powers, is associative::
1413 sage: set_random_seed()
1414 sage: x0 = random_eja().random_element()
1415 sage: A = x0.subalgebra_generated_by()
1416 sage: x,y,z = A.random_elements(3)
1417 sage: (x*y)*z == x*(y*z)
1420 Squaring in the subalgebra should work the same as in
1423 sage: set_random_seed()
1424 sage: x = random_eja().random_element()
1425 sage: A = x.subalgebra_generated_by()
1426 sage: A(x^2) == A(x)*A(x)
1429 By definition, the subalgebra generated by the zero element is
1430 the one-dimensional algebra generated by the identity
1431 element... unless the original algebra was trivial, in which
1432 case the subalgebra is trivial too::
1434 sage: set_random_seed()
1435 sage: A = random_eja().zero().subalgebra_generated_by()
1436 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1440 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1441 A
= self
.parent().subalgebra(powers
,
1446 A
.one
.set_cache(A(self
.parent().one()))
1450 def subalgebra_idempotent(self
):
1452 Find an idempotent in the associative subalgebra I generate
1453 using Proposition 2.3.5 in Baes.
1457 sage: from mjo.eja.eja_algebra import random_eja
1461 Ensure that we can find an idempotent in a non-trivial algebra
1462 where there are non-nilpotent elements, or that we get the dumb
1463 solution in the trivial algebra::
1465 sage: set_random_seed()
1466 sage: J = random_eja()
1467 sage: x = J.random_element()
1468 sage: while x.is_nilpotent() and not J.is_trivial():
1469 ....: x = J.random_element()
1470 sage: c = x.subalgebra_idempotent()
1475 if self
.parent().is_trivial():
1478 if self
.is_nilpotent():
1479 raise ValueError("this only works with non-nilpotent elements!")
1481 J
= self
.subalgebra_generated_by()
1484 # The image of the matrix of left-u^m-multiplication
1485 # will be minimal for some natural number s...
1487 minimal_dim
= J
.dimension()
1488 for i
in range(1, minimal_dim
):
1489 this_dim
= (u
**i
).operator().matrix().image().dimension()
1490 if this_dim
< minimal_dim
:
1491 minimal_dim
= this_dim
1494 # Now minimal_matrix should correspond to the smallest
1495 # non-zero subspace in Baes's (or really, Koecher's)
1498 # However, we need to restrict the matrix to work on the
1499 # subspace... or do we? Can't we just solve, knowing that
1500 # A(c) = u^(s+1) should have a solution in the big space,
1503 # Beware, solve_right() means that we're using COLUMN vectors.
1504 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1506 A
= u_next
.operator().matrix()
1507 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1509 # Now c is the idempotent we want, but it still lives in the subalgebra.
1510 return c
.superalgebra_element()
1515 Return my trace, the sum of my eigenvalues.
1517 In a trivial algebra, however you want to look at it, the trace is
1518 an empty sum for which we declare the result to be zero.
1522 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1529 sage: J = TrivialEJA()
1530 sage: J.zero().trace()
1534 sage: J = JordanSpinEJA(3)
1535 sage: x = sum(J.gens())
1541 sage: J = HadamardEJA(5)
1542 sage: J.one().trace()
1547 The trace of an element is a real number::
1549 sage: set_random_seed()
1550 sage: J = random_eja()
1551 sage: J.random_element().trace() in RLF
1554 The trace is linear::
1556 sage: set_random_seed()
1557 sage: J = random_eja()
1558 sage: x,y = J.random_elements(2)
1559 sage: alpha = J.base_ring().random_element()
1560 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1568 # Special case for the trivial algebra where
1569 # the trace is an empty sum.
1570 return P
.base_ring().zero()
1572 p
= P
._charpoly
_coefficients
()[r
-1]
1573 # The _charpoly_coeff function already adds the factor of
1574 # -1 to ensure that _charpoly_coeff(r-1) is really what
1575 # appears in front of t^{r-1} in the charpoly. However,
1576 # we want the negative of THAT for the trace.
1577 return -p(*self
.to_vector())
1580 def trace_inner_product(self
, other
):
1582 Return the trace inner product of myself and ``other``.
1586 sage: from mjo.eja.eja_algebra import random_eja
1590 The trace inner product is commutative, bilinear, and associative::
1592 sage: set_random_seed()
1593 sage: J = random_eja()
1594 sage: x,y,z = J.random_elements(3)
1596 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1599 sage: a = J.base_ring().random_element();
1600 sage: actual = (a*(x+z)).trace_inner_product(y)
1601 sage: expected = ( a*x.trace_inner_product(y) +
1602 ....: a*z.trace_inner_product(y) )
1603 sage: actual == expected
1605 sage: actual = x.trace_inner_product(a*(y+z))
1606 sage: expected = ( a*x.trace_inner_product(y) +
1607 ....: a*x.trace_inner_product(z) )
1608 sage: actual == expected
1611 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1615 if not other
in self
.parent():
1616 raise TypeError("'other' must live in the same algebra")
1618 return (self
*other
).trace()
1621 def trace_norm(self
):
1623 The norm of this element with respect to :meth:`trace_inner_product`.
1627 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1632 sage: J = HadamardEJA(2)
1633 sage: x = sum(J.gens())
1634 sage: x.trace_norm()
1639 sage: J = JordanSpinEJA(4)
1640 sage: x = sum(J.gens())
1641 sage: x.trace_norm()
1645 return self
.trace_inner_product(self
).sqrt()
1648 class CartesianProductEJAElement(FiniteDimensionalEJAElement
):
1651 Compute the determinant of this product-element using the
1652 determianants of its factors.
1654 This result Follows from the spectral decomposition of (say)
1655 the pair `(x,y)` in terms of the Jordan frame `\left\{ (c_1,
1656 0),(c_2, 0),...,(0,d_1),(0,d_2),... \right\}.
1658 from sage
.misc
.misc_c
import prod
1659 return prod( f
.det() for f
in self
.cartesian_factors() )
1661 def to_matrix(self
):
1662 # An override is necessary to call our custom _scale().
1663 B
= self
.parent().matrix_basis()
1664 W
= self
.parent().matrix_space()
1666 # Aaaaand linear combinations don't work in Cartesian
1667 # product spaces, even though they provide a method with
1668 # that name. This is hidden behind an "if" because the
1669 # _scale() function is slow.
1670 pairs
= zip(B
, self
.to_vector())
1671 return W
.sum( _scale(b
, alpha
) for (b
,alpha
) in pairs
)