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 (random_eja,
139 The rank of `R^3` is three, and the minimal polynomial of
140 the identity element is `(t-1)` from which it follows that
141 the characteristic polynomial should be `(t-1)^3`::
143 sage: J = HadamardEJA(3)
144 sage: J.one().characteristic_polynomial()
145 t^3 - 3*t^2 + 3*t - 1
147 Likewise, the characteristic of the zero element in the
148 rank-three algebra `R^{n}` should be `t^{3}`::
150 sage: J = HadamardEJA(3)
151 sage: J.zero().characteristic_polynomial()
156 The characteristic polynomial of an element should evaluate
157 to zero on that element::
159 sage: set_random_seed()
160 sage: x = random_eja().random_element()
161 sage: p = x.characteristic_polynomial()
162 sage: x.apply_univariate_polynomial(p).is_zero()
165 The characteristic polynomials of the zero and unit elements
166 should be what we think they are in a subalgebra, too::
168 sage: J = HadamardEJA(3)
169 sage: p1 = J.one().characteristic_polynomial()
170 sage: q1 = J.zero().characteristic_polynomial()
171 sage: b0,b1,b2 = J.gens()
172 sage: A = (b0 + 2*b1 + 3*b2).subalgebra_generated_by() # dim 3
173 sage: p2 = A.one().characteristic_polynomial()
174 sage: q2 = A.zero().characteristic_polynomial()
181 p
= self
.parent().characteristic_polynomial_of()
182 return p(*self
.to_vector())
185 def inner_product(self
, other
):
187 Return the parent algebra's inner product of myself and ``other``.
191 sage: from mjo.eja.eja_algebra import (
192 ....: ComplexHermitianEJA,
194 ....: QuaternionHermitianEJA,
195 ....: RealSymmetricEJA,
200 The inner product in the Jordan spin algebra is the usual
201 inner product on `R^n` (this example only works because the
202 basis for the Jordan algebra is the standard basis in `R^n`)::
204 sage: J = JordanSpinEJA(3)
205 sage: x = vector(QQ,[1,2,3])
206 sage: y = vector(QQ,[4,5,6])
207 sage: x.inner_product(y)
209 sage: J.from_vector(x).inner_product(J.from_vector(y))
212 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
213 multiplication is the usual matrix multiplication in `S^n`,
214 so the inner product of the identity matrix with itself
217 sage: J = RealSymmetricEJA(3)
218 sage: J.one().inner_product(J.one())
221 Likewise, the inner product on `C^n` is `<X,Y> =
222 Re(trace(X*Y))`, where we must necessarily take the real
223 part because the product of Hermitian matrices may not be
226 sage: J = ComplexHermitianEJA(3)
227 sage: J.one().inner_product(J.one())
230 Ditto for the quaternions::
232 sage: J = QuaternionHermitianEJA(2)
233 sage: J.one().inner_product(J.one())
238 Ensure that we can always compute an inner product, and that
239 it gives us back a real number::
241 sage: set_random_seed()
242 sage: J = random_eja()
243 sage: x,y = J.random_elements(2)
244 sage: x.inner_product(y) in RLF
250 raise TypeError("'other' must live in the same algebra")
252 return P
.inner_product(self
, other
)
255 def operator_commutes_with(self
, other
):
257 Return whether or not this element operator-commutes
262 sage: from mjo.eja.eja_algebra import random_eja
266 The definition of a Jordan algebra says that any element
267 operator-commutes with its square::
269 sage: set_random_seed()
270 sage: x = random_eja().random_element()
271 sage: x.operator_commutes_with(x^2)
276 Test Lemma 1 from Chapter III of Koecher::
278 sage: set_random_seed()
279 sage: u,v = random_eja().random_elements(2)
280 sage: lhs = u.operator_commutes_with(u*v)
281 sage: rhs = v.operator_commutes_with(u^2)
285 Test the first polarization identity from my notes, Koecher
286 Chapter III, or from Baes (2.3)::
288 sage: set_random_seed()
289 sage: x,y = random_eja().random_elements(2)
290 sage: Lx = x.operator()
291 sage: Ly = y.operator()
292 sage: Lxx = (x*x).operator()
293 sage: Lxy = (x*y).operator()
294 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
297 Test the second polarization identity from my notes or from
300 sage: set_random_seed()
301 sage: x,y,z = random_eja().random_elements(3)
302 sage: Lx = x.operator()
303 sage: Ly = y.operator()
304 sage: Lz = z.operator()
305 sage: Lzy = (z*y).operator()
306 sage: Lxy = (x*y).operator()
307 sage: Lxz = (x*z).operator()
308 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
311 Test the third polarization identity from my notes or from
314 sage: set_random_seed()
315 sage: u,y,z = random_eja().random_elements(3)
316 sage: Lu = u.operator()
317 sage: Ly = y.operator()
318 sage: Lz = z.operator()
319 sage: Lzy = (z*y).operator()
320 sage: Luy = (u*y).operator()
321 sage: Luz = (u*z).operator()
322 sage: Luyz = (u*(y*z)).operator()
323 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
324 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
325 sage: bool(lhs == rhs)
329 if not other
in self
.parent():
330 raise TypeError("'other' must live in the same algebra")
339 Return my determinant, the product of my eigenvalues.
343 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
345 ....: RealSymmetricEJA,
346 ....: ComplexHermitianEJA,
351 sage: J = JordanSpinEJA(2)
352 sage: x = sum( J.gens() )
358 sage: J = JordanSpinEJA(3)
359 sage: x = sum( J.gens() )
363 The determinant of the sole element in the rank-zero trivial
364 algebra is ``1``, by three paths of reasoning. First, its
365 characteristic polynomial is a constant ``1``, so the constant
366 term in that polynomial is ``1``. Second, the characteristic
367 polynomial evaluated at zero is again ``1``. And finally, the
368 (empty) product of its eigenvalues is likewise just unity::
370 sage: J = TrivialEJA()
376 An element is invertible if and only if its determinant is
379 sage: set_random_seed()
380 sage: x = random_eja().random_element()
381 sage: x.is_invertible() == (x.det() != 0)
384 Ensure that the determinant is multiplicative on an associative
385 subalgebra as in Faraut and Korányi's Proposition II.2.2::
387 sage: set_random_seed()
388 sage: J = random_eja().random_element().subalgebra_generated_by()
389 sage: x,y = J.random_elements(2)
390 sage: (x*y).det() == x.det()*y.det()
393 The determinant in real matrix algebras is the usual determinant::
395 sage: set_random_seed()
396 sage: X = matrix.random(QQ,3)
398 sage: J1 = RealSymmetricEJA(3)
399 sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
400 sage: expected = X.det()
401 sage: actual1 = J1(X).det()
402 sage: actual2 = J2(X).det()
403 sage: actual1 == expected
405 sage: actual2 == expected
413 # Special case, since we don't get the a0=1
414 # coefficient when the rank of the algebra
416 return P
.base_ring().one()
418 p
= P
._charpoly
_coefficients
()[0]
419 # The _charpoly_coeff function already adds the factor of -1
420 # to ensure that _charpoly_coefficients()[0] is really what
421 # appears in front of t^{0} in the charpoly. However, we want
422 # (-1)^r times THAT for the determinant.
423 return ((-1)**r
)*p(*self
.to_vector())
429 Return the Jordan-multiplicative inverse of this element.
433 In general we appeal to the quadratic representation as in
434 Koecher's Theorem 12 in Chapter III, Section 5. But if the
435 parent algebra's "characteristic polynomial of" coefficients
436 happen to be cached, then we use Proposition II.2.4 in Faraut
437 and Korányi which gives a formula for the inverse based on the
438 characteristic polynomial and the Cayley-Hamilton theorem for
439 Euclidean Jordan algebras::
443 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
449 The inverse in the spin factor algebra is given in Alizadeh's
452 sage: set_random_seed()
453 sage: J = JordanSpinEJA.random_instance()
454 sage: x = J.random_element()
455 sage: while not x.is_invertible():
456 ....: x = J.random_element()
457 sage: x_vec = x.to_vector()
459 sage: x_bar = x_vec[1:]
460 sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
461 sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
462 sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
463 sage: x.inverse() == J.from_vector(x_inverse)
466 Trying to invert a non-invertible element throws an error:
468 sage: JordanSpinEJA(3).zero().inverse()
469 Traceback (most recent call last):
471 ZeroDivisionError: element is not invertible
475 The identity element is its own inverse::
477 sage: set_random_seed()
478 sage: J = random_eja()
479 sage: J.one().inverse() == J.one()
482 If an element has an inverse, it acts like one::
484 sage: set_random_seed()
485 sage: J = random_eja()
486 sage: x = J.random_element()
487 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
490 The inverse of the inverse is what we started with::
492 sage: set_random_seed()
493 sage: J = random_eja()
494 sage: x = J.random_element()
495 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
498 Proposition II.2.3 in Faraut and Korányi says that the inverse
499 of an element is the inverse of its left-multiplication operator
500 applied to the algebra's identity, when that inverse exists::
502 sage: set_random_seed()
503 sage: J = random_eja()
504 sage: x = J.random_element()
505 sage: (not x.operator().is_invertible()) or (
506 ....: x.operator().inverse()(J.one()) == x.inverse() )
509 Check that the fast (cached) and slow algorithms give the same
512 sage: set_random_seed() # long time
513 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
514 sage: x = J.random_element() # long time
515 sage: while not x.is_invertible(): # long time
516 ....: x = J.random_element() # long time
517 sage: slow = x.inverse() # long time
518 sage: _ = J._charpoly_coefficients() # long time
519 sage: fast = x.inverse() # long time
520 sage: slow == fast # long time
523 not_invertible_msg
= "element is not invertible"
524 if self
.parent()._charpoly
_coefficients
.is_in_cache():
525 # We can invert using our charpoly if it will be fast to
526 # compute. If the coefficients are cached, our rank had
528 if self
.det().is_zero():
529 raise ZeroDivisionError(not_invertible_msg
)
530 r
= self
.parent().rank()
531 a
= self
.characteristic_polynomial().coefficients(sparse
=False)
532 return (-1)**(r
+1)*sum(a
[i
+1]*self
**i
for i
in range(r
))/self
.det()
535 inv
= (~self
.quadratic_representation())(self
)
536 self
.is_invertible
.set_cache(True)
538 except ZeroDivisionError:
539 self
.is_invertible
.set_cache(False)
540 raise ZeroDivisionError(not_invertible_msg
)
544 def is_invertible(self
):
546 Return whether or not this element is invertible.
550 If computing my determinant will be fast, we do so and compare
551 with zero (Proposition II.2.4 in Faraut and
552 Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
553 reduces the problem to the invertibility of my quadratic
558 sage: from mjo.eja.eja_algebra import random_eja
562 The identity element is always invertible::
564 sage: set_random_seed()
565 sage: J = random_eja()
566 sage: J.one().is_invertible()
569 The zero element is never invertible in a non-trivial algebra::
571 sage: set_random_seed()
572 sage: J = random_eja()
573 sage: (not J.is_trivial()) and J.zero().is_invertible()
576 Test that the fast (cached) and slow algorithms give the same
579 sage: set_random_seed() # long time
580 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
581 sage: x = J.random_element() # long time
582 sage: slow = x.is_invertible() # long time
583 sage: _ = J._charpoly_coefficients() # long time
584 sage: fast = x.is_invertible() # long time
585 sage: slow == fast # long time
589 if self
.parent().is_trivial():
594 if self
.parent()._charpoly
_coefficients
.is_in_cache():
595 # The determinant will be quicker than inverting the
596 # quadratic representation, most likely.
597 return (not self
.det().is_zero())
599 # The easiest way to determine if I'm invertible is to try.
601 inv
= (~self
.quadratic_representation())(self
)
602 self
.inverse
.set_cache(inv
)
604 except ZeroDivisionError:
608 def is_primitive_idempotent(self
):
610 Return whether or not this element is a primitive (or minimal)
613 A primitive idempotent is a non-zero idempotent that is not
614 the sum of two other non-zero idempotents. Remark 2.7.15 in
615 Baes shows that this is what he refers to as a "minimal
618 An element of a Euclidean Jordan algebra is a minimal idempotent
619 if it :meth:`is_idempotent` and if its Peirce subalgebra
620 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
625 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
626 ....: RealSymmetricEJA,
632 This method is sloooooow.
636 The spectral decomposition of a non-regular element should always
637 contain at least one non-minimal idempotent::
639 sage: J = RealSymmetricEJA(3)
640 sage: x = sum(J.gens())
643 sage: [ c.is_primitive_idempotent()
644 ....: for (l,c) in x.spectral_decomposition() ]
647 On the other hand, the spectral decomposition of a regular
648 element should always be in terms of minimal idempotents::
650 sage: J = JordanSpinEJA(4)
651 sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
654 sage: [ c.is_primitive_idempotent()
655 ....: for (l,c) in x.spectral_decomposition() ]
660 The identity element is minimal only in an EJA of rank one::
662 sage: set_random_seed()
663 sage: J = random_eja()
664 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
667 A non-idempotent cannot be a minimal idempotent::
669 sage: set_random_seed()
670 sage: J = JordanSpinEJA(4)
671 sage: x = J.random_element()
672 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
675 Proposition 2.7.19 in Baes says that an element is a minimal
676 idempotent if and only if it's idempotent with trace equal to
679 sage: set_random_seed()
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: set_random_seed()
690 sage: J = random_eja()
691 sage: J.zero().is_idempotent()
693 sage: J.zero().is_primitive_idempotent()
696 As a consequence of the fact that primitive idempotents must
697 be non-zero, there are no primitive idempotents in a trivial
698 Euclidean Jordan algebra::
700 sage: J = TrivialEJA()
701 sage: J.one().is_idempotent()
703 sage: J.one().is_primitive_idempotent()
707 if not self
.is_idempotent():
713 (_
,_
,J1
) = self
.parent().peirce_decomposition(self
)
714 return (J1
.dimension() == 1)
717 def is_nilpotent(self
):
719 Return whether or not some power of this element is zero.
723 We use Theorem 5 in Chapter III of Koecher, which says that
724 an element ``x`` is nilpotent if and only if ``x.operator()``
725 is nilpotent. And it is a basic fact of linear algebra that
726 an operator on an `n`-dimensional space is nilpotent if and
727 only if, when raised to the `n`th power, it equals the zero
728 operator (for example, see Axler Corollary 8.8).
732 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
737 sage: J = JordanSpinEJA(3)
738 sage: x = sum(J.gens())
739 sage: x.is_nilpotent()
744 The identity element is never nilpotent, except in a trivial EJA::
746 sage: set_random_seed()
747 sage: J = random_eja()
748 sage: J.one().is_nilpotent() and not J.is_trivial()
751 The additive identity is always nilpotent::
753 sage: set_random_seed()
754 sage: random_eja().zero().is_nilpotent()
759 zero_operator
= P
.zero().operator()
760 return self
.operator()**P
.dimension() == zero_operator
763 def is_regular(self
):
765 Return whether or not this is a regular element.
769 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
774 The identity element always has degree one, but any element
775 linearly-independent from it is regular::
777 sage: J = JordanSpinEJA(5)
778 sage: J.one().is_regular()
780 sage: b0, b1, b2, b3, b4 = J.gens()
783 sage: for x in J.gens():
784 ....: (J.one() + x).is_regular()
793 The zero element should never be regular, unless the parent
794 algebra has dimension less than or equal to one::
796 sage: set_random_seed()
797 sage: J = random_eja()
798 sage: J.dimension() <= 1 or not J.zero().is_regular()
801 The unit element isn't regular unless the algebra happens to
802 consist of only its scalar multiples::
804 sage: set_random_seed()
805 sage: J = random_eja()
806 sage: J.dimension() <= 1 or not J.one().is_regular()
810 return self
.degree() == self
.parent().rank()
815 Return the degree of this element, which is defined to be
816 the degree of its minimal polynomial.
824 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
829 sage: J = JordanSpinEJA(4)
830 sage: J.one().degree()
832 sage: b0,b1,b2,b3 = J.gens()
833 sage: (b0 - b1).degree()
836 In the spin factor algebra (of rank two), all elements that
837 aren't multiples of the identity are regular::
839 sage: set_random_seed()
840 sage: J = JordanSpinEJA.random_instance()
841 sage: n = J.dimension()
842 sage: x = J.random_element()
843 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
848 The zero and unit elements are both of degree one in nontrivial
851 sage: set_random_seed()
852 sage: J = random_eja()
853 sage: d = J.zero().degree()
854 sage: (J.is_trivial() and d == 0) or d == 1
856 sage: d = J.one().degree()
857 sage: (J.is_trivial() and d == 0) or d == 1
860 Our implementation agrees with the definition::
862 sage: set_random_seed()
863 sage: x = random_eja().random_element()
864 sage: x.degree() == x.minimal_polynomial().degree()
868 n
= self
.parent().dimension()
871 # The minimal polynomial is an empty product, i.e. the
872 # constant polynomial "1" having degree zero.
875 # The minimal polynomial of zero in a nontrivial algebra
876 # is "t", and is of degree one.
879 # If this is a nonzero element of a nontrivial algebra, it
880 # has degree at least one. It follows that, in an algebra
881 # of dimension one, the degree must be actually one.
884 # BEWARE: The subalgebra_generated_by() method uses the result
885 # of this method to construct a basis for the subalgebra. That
886 # means, in particular, that we cannot implement this method
887 # as ``self.subalgebra_generated_by().dimension()``.
889 # Algorithm: keep appending (vector representations of) powers
890 # self as rows to a matrix and echelonizing it. When its rank
891 # stops increasing, we've reached a redundancy.
893 # Given the special cases above, we can assume that "self" is
894 # nonzero, the algebra is nontrivial, and that its dimension
896 M
= matrix([(self
.parent().one()).to_vector()])
899 # Specifying the row-reduction algorithm can e.g. help over
900 # AA because it avoids the RecursionError that gets thrown
901 # when we have to look too hard for a root.
903 # Beware: QQ supports an entirely different set of "algorithm"
904 # keywords than do AA and RR.
906 from sage
.rings
.all
import QQ
907 if self
.parent().base_ring() is not QQ
:
908 algo
= "scaled_partial_pivoting"
911 M
= matrix(M
.rows() + [(self
**d
).to_vector()])
914 if new_rank
== old_rank
:
923 def left_matrix(self
):
925 Our parent class defines ``left_matrix`` and ``matrix``
926 methods whose names are misleading. We don't want them.
928 raise NotImplementedError("use operator().matrix() instead")
933 def minimal_polynomial(self
):
935 Return the minimal polynomial of this element,
936 as a function of the variable `t`.
940 We restrict ourselves to the associative subalgebra
941 generated by this element, and then return the minimal
942 polynomial of this element's operator matrix (in that
943 subalgebra). This works by Baes Proposition 2.3.16.
947 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
948 ....: RealSymmetricEJA,
954 Keeping in mind that the polynomial ``1`` evaluates the identity
955 element (also the zero element) of the trivial algebra, it is clear
956 that the polynomial ``1`` is the minimal polynomial of the only
957 element in a trivial algebra::
959 sage: J = TrivialEJA()
960 sage: J.one().minimal_polynomial()
962 sage: J.zero().minimal_polynomial()
967 The minimal polynomial of the identity and zero elements are
968 always the same, except in trivial algebras where the minimal
969 polynomial of the unit/zero element is ``1``::
971 sage: set_random_seed()
972 sage: J = random_eja()
973 sage: mu = J.one().minimal_polynomial()
974 sage: t = mu.parent().gen()
975 sage: mu + int(J.is_trivial())*(t-2)
977 sage: mu = J.zero().minimal_polynomial()
978 sage: t = mu.parent().gen()
979 sage: mu + int(J.is_trivial())*(t-1)
982 The degree of an element is (by one definition) the degree
983 of its minimal polynomial::
985 sage: set_random_seed()
986 sage: x = random_eja().random_element()
987 sage: x.degree() == x.minimal_polynomial().degree()
990 The minimal polynomial and the characteristic polynomial coincide
991 and are known (see Alizadeh, Example 11.11) for all elements of
992 the spin factor algebra that aren't scalar multiples of the
993 identity. We require the dimension of the algebra to be at least
994 two here so that said elements actually exist::
996 sage: set_random_seed()
997 sage: d_max = JordanSpinEJA._max_random_instance_dimension()
998 sage: n = ZZ.random_element(2, max(2,d_max))
999 sage: J = JordanSpinEJA(n)
1000 sage: y = J.random_element()
1001 sage: while y == y.coefficient(0)*J.one():
1002 ....: y = J.random_element()
1003 sage: y0 = y.to_vector()[0]
1004 sage: y_bar = y.to_vector()[1:]
1005 sage: actual = y.minimal_polynomial()
1006 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1007 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1008 sage: bool(actual == expected)
1011 The minimal polynomial should always kill its element::
1013 sage: set_random_seed()
1014 sage: x = random_eja().random_element()
1015 sage: p = x.minimal_polynomial()
1016 sage: x.apply_univariate_polynomial(p)
1019 The minimal polynomial is invariant under a change of basis,
1020 and in particular, a re-scaling of the basis::
1022 sage: set_random_seed()
1023 sage: d_max = RealSymmetricEJA._max_random_instance_dimension()
1024 sage: n = ZZ.random_element(1, d_max)
1025 sage: J1 = RealSymmetricEJA(n)
1026 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
1027 sage: X = random_matrix(AA,n)
1028 sage: X = X*X.transpose()
1031 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
1036 # Pretty sure we know what the minimal polynomial of
1037 # the zero operator is going to be. This ensures
1038 # consistency of e.g. the polynomial variable returned
1039 # in the "normal" case without us having to think about it.
1040 return self
.operator().minimal_polynomial()
1042 # If we don't orthonormalize the subalgebra's basis, then the
1043 # first two monomials in the subalgebra will be self^0 and
1044 # self^1... assuming that self^1 is not a scalar multiple of
1045 # self^0 (the unit element). We special case these to avoid
1046 # having to solve a system to coerce self into the subalgebra.
1047 A
= self
.subalgebra_generated_by(orthonormalize
=False)
1049 if A
.dimension() == 1:
1050 # Does a solve to find the scalar multiple alpha such that
1051 # alpha*unit = self. We have to do this because the basis
1052 # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
1053 unit
= self
.parent().one()
1054 alpha
= self
.to_vector() / unit
.to_vector()
1055 return (unit
.operator()*alpha
).minimal_polynomial()
1057 # If the dimension of the subalgebra is >= 2, then we just
1058 # use the second basis element.
1059 return A
.monomial(1).operator().minimal_polynomial()
1063 def to_matrix(self
):
1065 Return an (often more natural) representation of this element as a
1068 Every finite-dimensional Euclidean Jordan Algebra is a direct
1069 sum of five simple algebras, four of which comprise Hermitian
1070 matrices. This method returns a "natural" matrix
1071 representation of this element as either a Hermitian matrix or
1076 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1078 ....: QuaternionHermitianEJA,
1079 ....: RealSymmetricEJA)
1083 sage: J = ComplexHermitianEJA(3)
1086 sage: J.one().to_matrix()
1097 sage: J = QuaternionHermitianEJA(2)
1100 sage: J.one().to_matrix()
1107 This also works in Cartesian product algebras::
1109 sage: J1 = HadamardEJA(1)
1110 sage: J2 = RealSymmetricEJA(2)
1111 sage: J = cartesian_product([J1,J2])
1112 sage: x = sum(J.gens())
1113 sage: x.to_matrix()[0]
1115 sage: x.to_matrix()[1]
1116 [ 1 0.7071067811865475?]
1117 [0.7071067811865475? 1]
1120 B
= self
.parent().matrix_basis()
1121 W
= self
.parent().matrix_space()
1123 if hasattr(W
, 'cartesian_factors'):
1124 # Aaaaand linear combinations don't work in Cartesian
1125 # product spaces, even though they provide a method with
1126 # that name. This is hidden behind an "if" because the
1127 # _scale() function is slow.
1128 pairs
= zip(B
, self
.to_vector())
1129 return W
.sum( _scale(b
, alpha
) for (b
,alpha
) in pairs
)
1131 # This is just a manual "from_vector()", but of course
1132 # matrix spaces aren't vector spaces in sage, so they
1133 # don't have a from_vector() method.
1134 return W
.linear_combination( zip(B
, self
.to_vector()) )
1140 The norm of this element with respect to :meth:`inner_product`.
1144 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1149 sage: J = HadamardEJA(2)
1150 sage: x = sum(J.gens())
1156 sage: J = JordanSpinEJA(4)
1157 sage: x = sum(J.gens())
1162 return self
.inner_product(self
).sqrt()
1167 Return the left-multiplication-by-this-element
1168 operator on the ambient algebra.
1172 sage: from mjo.eja.eja_algebra import random_eja
1176 sage: set_random_seed()
1177 sage: J = random_eja()
1178 sage: x,y = J.random_elements(2)
1179 sage: x.operator()(y) == x*y
1181 sage: y.operator()(x) == x*y
1186 left_mult_by_self
= lambda y
: self
*y
1187 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1188 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1191 def quadratic_representation(self
, other
=None):
1193 Return the quadratic representation of this element.
1197 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1202 The explicit form in the spin factor algebra is given by
1203 Alizadeh's Example 11.12::
1205 sage: set_random_seed()
1206 sage: x = JordanSpinEJA.random_instance().random_element()
1207 sage: x_vec = x.to_vector()
1208 sage: Q = matrix.identity(x.base_ring(), 0)
1209 sage: n = x_vec.degree()
1212 ....: x_bar = x_vec[1:]
1213 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1214 ....: B = 2*x0*x_bar.row()
1215 ....: C = 2*x0*x_bar.column()
1216 ....: D = matrix.identity(x.base_ring(), n-1)
1217 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1218 ....: D = D + 2*x_bar.tensor_product(x_bar)
1219 ....: Q = matrix.block(2,2,[A,B,C,D])
1220 sage: Q == x.quadratic_representation().matrix()
1223 Test all of the properties from Theorem 11.2 in Alizadeh::
1225 sage: set_random_seed()
1226 sage: J = random_eja()
1227 sage: x,y = J.random_elements(2)
1228 sage: Lx = x.operator()
1229 sage: Lxx = (x*x).operator()
1230 sage: Qx = x.quadratic_representation()
1231 sage: Qy = y.quadratic_representation()
1232 sage: Qxy = x.quadratic_representation(y)
1233 sage: Qex = J.one().quadratic_representation(x)
1234 sage: n = ZZ.random_element(10)
1235 sage: Qxn = (x^n).quadratic_representation()
1239 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1242 Property 2 (multiply on the right for :trac:`28272`):
1244 sage: alpha = J.base_ring().random_element()
1245 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1250 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1253 sage: not x.is_invertible() or (
1256 ....: x.inverse().quadratic_representation() )
1259 sage: Qxy(J.one()) == x*y
1264 sage: not x.is_invertible() or (
1265 ....: x.quadratic_representation(x.inverse())*Qx
1266 ....: == Qx*x.quadratic_representation(x.inverse()) )
1269 sage: not x.is_invertible() or (
1270 ....: x.quadratic_representation(x.inverse())*Qx
1272 ....: 2*Lx*Qex - Qx )
1275 sage: 2*Lx*Qex - Qx == Lxx
1280 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1290 sage: not x.is_invertible() or (
1291 ....: Qx*x.inverse().operator() == Lx )
1296 sage: not x.operator_commutes_with(y) or (
1297 ....: Qx(y)^n == Qxn(y^n) )
1303 elif not other
in self
.parent():
1304 raise TypeError("'other' must live in the same algebra")
1307 M
= other
.operator()
1308 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1312 def spectral_decomposition(self
):
1314 Return the unique spectral decomposition of this element.
1318 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1319 element's left-multiplication-by operator to the subalgebra it
1320 generates. We then compute the spectral decomposition of that
1321 operator, and the spectral projectors we get back must be the
1322 left-multiplication-by operators for the idempotents we
1323 seek. Thus applying them to the identity element gives us those
1326 Since the eigenvalues are required to be distinct, we take
1327 the spectral decomposition of the zero element to be zero
1328 times the identity element of the algebra (which is idempotent,
1333 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1337 The spectral decomposition of the identity is ``1`` times itself,
1338 and the spectral decomposition of zero is ``0`` times the identity::
1340 sage: J = RealSymmetricEJA(3)
1343 sage: J.one().spectral_decomposition()
1345 sage: J.zero().spectral_decomposition()
1350 sage: J = RealSymmetricEJA(4)
1351 sage: x = sum(J.gens())
1352 sage: sd = x.spectral_decomposition()
1357 sage: c0.inner_product(c1) == 0
1359 sage: c0.is_idempotent()
1361 sage: c1.is_idempotent()
1363 sage: c0 + c1 == J.one()
1365 sage: l0*c0 + l1*c1 == x
1368 The spectral decomposition should work in subalgebras, too::
1370 sage: J = RealSymmetricEJA(4)
1371 sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
1372 sage: A = 2*b5 - 2*b8
1373 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1374 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1375 sage: (f0, f1, f2) = J1.gens()
1376 sage: f0.spectral_decomposition()
1380 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1382 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1383 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1386 def subalgebra_generated_by(self
, **kwargs
):
1388 Return the associative subalgebra of the parent EJA generated
1391 Since our parent algebra is unital, we want "subalgebra" to mean
1392 "unital subalgebra" as well; thus the subalgebra that an element
1393 generates will itself be a Euclidean Jordan algebra after
1394 restricting the algebra operations appropriately. This is the
1395 subalgebra that Faraut and Korányi work with in section II.2, for
1400 sage: from mjo.eja.eja_algebra import (random_eja,
1402 ....: RealSymmetricEJA)
1406 We can create subalgebras of Cartesian product EJAs that are not
1407 themselves Cartesian product EJAs (they're just "regular" EJAs)::
1409 sage: J1 = HadamardEJA(3)
1410 sage: J2 = RealSymmetricEJA(2)
1411 sage: J = cartesian_product([J1,J2])
1412 sage: J.one().subalgebra_generated_by()
1413 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
1417 This subalgebra, being composed of only powers, is associative::
1419 sage: set_random_seed()
1420 sage: x0 = random_eja().random_element()
1421 sage: A = x0.subalgebra_generated_by()
1422 sage: x,y,z = A.random_elements(3)
1423 sage: (x*y)*z == x*(y*z)
1426 Squaring in the subalgebra should work the same as in
1429 sage: set_random_seed()
1430 sage: x = random_eja().random_element()
1431 sage: A = x.subalgebra_generated_by()
1432 sage: A(x^2) == A(x)*A(x)
1435 By definition, the subalgebra generated by the zero element is
1436 the one-dimensional algebra generated by the identity
1437 element... unless the original algebra was trivial, in which
1438 case the subalgebra is trivial too::
1440 sage: set_random_seed()
1441 sage: A = random_eja().zero().subalgebra_generated_by()
1442 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1446 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1447 A
= self
.parent().subalgebra(powers
,
1452 A
.one
.set_cache(A(self
.parent().one()))
1456 def subalgebra_idempotent(self
):
1458 Find an idempotent in the associative subalgebra I generate
1459 using Proposition 2.3.5 in Baes.
1463 sage: from mjo.eja.eja_algebra import random_eja
1467 Ensure that we can find an idempotent in a non-trivial algebra
1468 where there are non-nilpotent elements, or that we get the dumb
1469 solution in the trivial algebra::
1471 sage: set_random_seed()
1472 sage: J = random_eja()
1473 sage: x = J.random_element()
1474 sage: while x.is_nilpotent() and not J.is_trivial():
1475 ....: x = J.random_element()
1476 sage: c = x.subalgebra_idempotent()
1481 if self
.parent().is_trivial():
1484 if self
.is_nilpotent():
1485 raise ValueError("this only works with non-nilpotent elements!")
1487 J
= self
.subalgebra_generated_by()
1490 # The image of the matrix of left-u^m-multiplication
1491 # will be minimal for some natural number s...
1493 minimal_dim
= J
.dimension()
1494 for i
in range(1, minimal_dim
):
1495 this_dim
= (u
**i
).operator().matrix().image().dimension()
1496 if this_dim
< minimal_dim
:
1497 minimal_dim
= this_dim
1500 # Now minimal_matrix should correspond to the smallest
1501 # non-zero subspace in Baes's (or really, Koecher's)
1504 # However, we need to restrict the matrix to work on the
1505 # subspace... or do we? Can't we just solve, knowing that
1506 # A(c) = u^(s+1) should have a solution in the big space,
1509 # Beware, solve_right() means that we're using COLUMN vectors.
1510 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1512 A
= u_next
.operator().matrix()
1513 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1515 # Now c is the idempotent we want, but it still lives in the subalgebra.
1516 return c
.superalgebra_element()
1521 Return my trace, the sum of my eigenvalues.
1523 In a trivial algebra, however you want to look at it, the trace is
1524 an empty sum for which we declare the result to be zero.
1528 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1535 sage: J = TrivialEJA()
1536 sage: J.zero().trace()
1540 sage: J = JordanSpinEJA(3)
1541 sage: x = sum(J.gens())
1547 sage: J = HadamardEJA(5)
1548 sage: J.one().trace()
1553 The trace of an element is a real number::
1555 sage: set_random_seed()
1556 sage: J = random_eja()
1557 sage: J.random_element().trace() in RLF
1560 The trace is linear::
1562 sage: set_random_seed()
1563 sage: J = random_eja()
1564 sage: x,y = J.random_elements(2)
1565 sage: alpha = J.base_ring().random_element()
1566 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1574 # Special case for the trivial algebra where
1575 # the trace is an empty sum.
1576 return P
.base_ring().zero()
1578 p
= P
._charpoly
_coefficients
()[r
-1]
1579 # The _charpoly_coeff function already adds the factor of
1580 # -1 to ensure that _charpoly_coeff(r-1) is really what
1581 # appears in front of t^{r-1} in the charpoly. However,
1582 # we want the negative of THAT for the trace.
1583 return -p(*self
.to_vector())
1586 def trace_inner_product(self
, other
):
1588 Return the trace inner product of myself and ``other``.
1592 sage: from mjo.eja.eja_algebra import random_eja
1596 The trace inner product is commutative, bilinear, and associative::
1598 sage: set_random_seed()
1599 sage: J = random_eja()
1600 sage: x,y,z = J.random_elements(3)
1602 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1605 sage: a = J.base_ring().random_element();
1606 sage: actual = (a*(x+z)).trace_inner_product(y)
1607 sage: expected = ( a*x.trace_inner_product(y) +
1608 ....: a*z.trace_inner_product(y) )
1609 sage: actual == expected
1611 sage: actual = x.trace_inner_product(a*(y+z))
1612 sage: expected = ( a*x.trace_inner_product(y) +
1613 ....: a*x.trace_inner_product(z) )
1614 sage: actual == expected
1617 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1621 if not other
in self
.parent():
1622 raise TypeError("'other' must live in the same algebra")
1624 return (self
*other
).trace()
1627 def trace_norm(self
):
1629 The norm of this element with respect to :meth:`trace_inner_product`.
1633 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1638 sage: J = HadamardEJA(2)
1639 sage: x = sum(J.gens())
1640 sage: x.trace_norm()
1645 sage: J = JordanSpinEJA(4)
1646 sage: x = sum(J.gens())
1647 sage: x.trace_norm()
1651 return self
.trace_inner_product(self
).sqrt()