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
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: d = ZZ.random_element(1, d_max)
1025 sage: n = RealSymmetricEJA._max_random_instance_size(d)
1026 sage: J1 = RealSymmetricEJA(n)
1027 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
1028 sage: X = random_matrix(AA,n)
1029 sage: X = X*X.transpose()
1032 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
1037 # Pretty sure we know what the minimal polynomial of
1038 # the zero operator is going to be. This ensures
1039 # consistency of e.g. the polynomial variable returned
1040 # in the "normal" case without us having to think about it.
1041 return self
.operator().minimal_polynomial()
1043 # If we don't orthonormalize the subalgebra's basis, then the
1044 # first two monomials in the subalgebra will be self^0 and
1045 # self^1... assuming that self^1 is not a scalar multiple of
1046 # self^0 (the unit element). We special case these to avoid
1047 # having to solve a system to coerce self into the subalgebra.
1048 A
= self
.subalgebra_generated_by(orthonormalize
=False)
1050 if A
.dimension() == 1:
1051 # Does a solve to find the scalar multiple alpha such that
1052 # alpha*unit = self. We have to do this because the basis
1053 # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
1054 unit
= self
.parent().one()
1055 alpha
= self
.to_vector() / unit
.to_vector()
1056 return (unit
.operator()*alpha
).minimal_polynomial()
1058 # If the dimension of the subalgebra is >= 2, then we just
1059 # use the second basis element.
1060 return A
.monomial(1).operator().minimal_polynomial()
1064 def to_matrix(self
):
1066 Return an (often more natural) representation of this element as a
1069 Every finite-dimensional Euclidean Jordan Algebra is a direct
1070 sum of five simple algebras, four of which comprise Hermitian
1071 matrices. This method returns a "natural" matrix
1072 representation of this element as either a Hermitian matrix or
1077 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1079 ....: QuaternionHermitianEJA,
1080 ....: RealSymmetricEJA)
1084 sage: J = ComplexHermitianEJA(3)
1087 sage: J.one().to_matrix()
1098 sage: J = QuaternionHermitianEJA(2)
1101 sage: J.one().to_matrix()
1108 This also works in Cartesian product algebras::
1110 sage: J1 = HadamardEJA(1)
1111 sage: J2 = RealSymmetricEJA(2)
1112 sage: J = cartesian_product([J1,J2])
1113 sage: x = sum(J.gens())
1114 sage: x.to_matrix()[0]
1116 sage: x.to_matrix()[1]
1117 [ 1 0.7071067811865475?]
1118 [0.7071067811865475? 1]
1121 B
= self
.parent().matrix_basis()
1122 W
= self
.parent().matrix_space()
1124 if hasattr(W
, 'cartesian_factors'):
1125 # Aaaaand linear combinations don't work in Cartesian
1126 # product spaces, even though they provide a method with
1127 # that name. This is hidden behind an "if" because the
1128 # _scale() function is slow.
1129 pairs
= zip(B
, self
.to_vector())
1130 return W
.sum( _scale(b
, alpha
) for (b
,alpha
) in pairs
)
1132 # This is just a manual "from_vector()", but of course
1133 # matrix spaces aren't vector spaces in sage, so they
1134 # don't have a from_vector() method.
1135 return W
.linear_combination( zip(B
, self
.to_vector()) )
1141 The norm of this element with respect to :meth:`inner_product`.
1145 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1150 sage: J = HadamardEJA(2)
1151 sage: x = sum(J.gens())
1157 sage: J = JordanSpinEJA(4)
1158 sage: x = sum(J.gens())
1163 return self
.inner_product(self
).sqrt()
1168 Return the left-multiplication-by-this-element
1169 operator on the ambient algebra.
1173 sage: from mjo.eja.eja_algebra import random_eja
1177 sage: set_random_seed()
1178 sage: J = random_eja()
1179 sage: x,y = J.random_elements(2)
1180 sage: x.operator()(y) == x*y
1182 sage: y.operator()(x) == x*y
1187 left_mult_by_self
= lambda y
: self
*y
1188 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1189 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1192 def quadratic_representation(self
, other
=None):
1194 Return the quadratic representation of this element.
1198 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1203 The explicit form in the spin factor algebra is given by
1204 Alizadeh's Example 11.12::
1206 sage: set_random_seed()
1207 sage: x = JordanSpinEJA.random_instance().random_element()
1208 sage: x_vec = x.to_vector()
1209 sage: Q = matrix.identity(x.base_ring(), 0)
1210 sage: n = x_vec.degree()
1213 ....: x_bar = x_vec[1:]
1214 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1215 ....: B = 2*x0*x_bar.row()
1216 ....: C = 2*x0*x_bar.column()
1217 ....: D = matrix.identity(x.base_ring(), n-1)
1218 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1219 ....: D = D + 2*x_bar.tensor_product(x_bar)
1220 ....: Q = matrix.block(2,2,[A,B,C,D])
1221 sage: Q == x.quadratic_representation().matrix()
1224 Test all of the properties from Theorem 11.2 in Alizadeh::
1226 sage: set_random_seed()
1227 sage: J = random_eja()
1228 sage: x,y = J.random_elements(2)
1229 sage: Lx = x.operator()
1230 sage: Lxx = (x*x).operator()
1231 sage: Qx = x.quadratic_representation()
1232 sage: Qy = y.quadratic_representation()
1233 sage: Qxy = x.quadratic_representation(y)
1234 sage: Qex = J.one().quadratic_representation(x)
1235 sage: n = ZZ.random_element(10)
1236 sage: Qxn = (x^n).quadratic_representation()
1240 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1243 Property 2 (multiply on the right for :trac:`28272`):
1245 sage: alpha = J.base_ring().random_element()
1246 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1251 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1254 sage: not x.is_invertible() or (
1257 ....: x.inverse().quadratic_representation() )
1260 sage: Qxy(J.one()) == x*y
1265 sage: not x.is_invertible() or (
1266 ....: x.quadratic_representation(x.inverse())*Qx
1267 ....: == Qx*x.quadratic_representation(x.inverse()) )
1270 sage: not x.is_invertible() or (
1271 ....: x.quadratic_representation(x.inverse())*Qx
1273 ....: 2*Lx*Qex - Qx )
1276 sage: 2*Lx*Qex - Qx == Lxx
1281 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1291 sage: not x.is_invertible() or (
1292 ....: Qx*x.inverse().operator() == Lx )
1297 sage: not x.operator_commutes_with(y) or (
1298 ....: Qx(y)^n == Qxn(y^n) )
1304 elif not other
in self
.parent():
1305 raise TypeError("'other' must live in the same algebra")
1308 M
= other
.operator()
1309 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1313 def spectral_decomposition(self
):
1315 Return the unique spectral decomposition of this element.
1319 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1320 element's left-multiplication-by operator to the subalgebra it
1321 generates. We then compute the spectral decomposition of that
1322 operator, and the spectral projectors we get back must be the
1323 left-multiplication-by operators for the idempotents we
1324 seek. Thus applying them to the identity element gives us those
1327 Since the eigenvalues are required to be distinct, we take
1328 the spectral decomposition of the zero element to be zero
1329 times the identity element of the algebra (which is idempotent,
1334 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1338 The spectral decomposition of the identity is ``1`` times itself,
1339 and the spectral decomposition of zero is ``0`` times the identity::
1341 sage: J = RealSymmetricEJA(3)
1344 sage: J.one().spectral_decomposition()
1346 sage: J.zero().spectral_decomposition()
1351 sage: J = RealSymmetricEJA(4)
1352 sage: x = sum(J.gens())
1353 sage: sd = x.spectral_decomposition()
1358 sage: c0.inner_product(c1) == 0
1360 sage: c0.is_idempotent()
1362 sage: c1.is_idempotent()
1364 sage: c0 + c1 == J.one()
1366 sage: l0*c0 + l1*c1 == x
1369 The spectral decomposition should work in subalgebras, too::
1371 sage: J = RealSymmetricEJA(4)
1372 sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
1373 sage: A = 2*b5 - 2*b8
1374 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1375 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1376 sage: (f0, f1, f2) = J1.gens()
1377 sage: f0.spectral_decomposition()
1378 [(0, 1.000000000000000?*c2), (1, 1.000000000000000?*c0)]
1381 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1383 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1384 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1387 def subalgebra_generated_by(self
, **kwargs
):
1389 Return the associative subalgebra of the parent EJA generated
1392 Since our parent algebra is unital, we want "subalgebra" to mean
1393 "unital subalgebra" as well; thus the subalgebra that an element
1394 generates will itself be a Euclidean Jordan algebra after
1395 restricting the algebra operations appropriately. This is the
1396 subalgebra that Faraut and Korányi work with in section II.2, for
1401 sage: from mjo.eja.eja_algebra import (random_eja,
1403 ....: RealSymmetricEJA)
1407 We can create subalgebras of Cartesian product EJAs that are not
1408 themselves Cartesian product EJAs (they're just "regular" EJAs)::
1410 sage: J1 = HadamardEJA(3)
1411 sage: J2 = RealSymmetricEJA(2)
1412 sage: J = cartesian_product([J1,J2])
1413 sage: J.one().subalgebra_generated_by()
1414 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
1418 This subalgebra, being composed of only powers, is associative::
1420 sage: set_random_seed()
1421 sage: x0 = random_eja().random_element()
1422 sage: A = x0.subalgebra_generated_by()
1423 sage: x,y,z = A.random_elements(3)
1424 sage: (x*y)*z == x*(y*z)
1427 Squaring in the subalgebra should work the same as in
1430 sage: set_random_seed()
1431 sage: x = random_eja().random_element()
1432 sage: A = x.subalgebra_generated_by()
1433 sage: A(x^2) == A(x)*A(x)
1436 By definition, the subalgebra generated by the zero element is
1437 the one-dimensional algebra generated by the identity
1438 element... unless the original algebra was trivial, in which
1439 case the subalgebra is trivial too::
1441 sage: set_random_seed()
1442 sage: A = random_eja().zero().subalgebra_generated_by()
1443 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1447 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1448 A
= self
.parent().subalgebra(powers
,
1453 A
.one
.set_cache(A(self
.parent().one()))
1457 def subalgebra_idempotent(self
):
1459 Find an idempotent in the associative subalgebra I generate
1460 using Proposition 2.3.5 in Baes.
1464 sage: from mjo.eja.eja_algebra import random_eja
1468 Ensure that we can find an idempotent in a non-trivial algebra
1469 where there are non-nilpotent elements, or that we get the dumb
1470 solution in the trivial algebra::
1472 sage: set_random_seed()
1473 sage: J = random_eja()
1474 sage: x = J.random_element()
1475 sage: while x.is_nilpotent() and not J.is_trivial():
1476 ....: x = J.random_element()
1477 sage: c = x.subalgebra_idempotent()
1482 if self
.parent().is_trivial():
1485 if self
.is_nilpotent():
1486 raise ValueError("this only works with non-nilpotent elements!")
1488 J
= self
.subalgebra_generated_by()
1491 # The image of the matrix of left-u^m-multiplication
1492 # will be minimal for some natural number s...
1494 minimal_dim
= J
.dimension()
1495 for i
in range(1, minimal_dim
):
1496 this_dim
= (u
**i
).operator().matrix().image().dimension()
1497 if this_dim
< minimal_dim
:
1498 minimal_dim
= this_dim
1501 # Now minimal_matrix should correspond to the smallest
1502 # non-zero subspace in Baes's (or really, Koecher's)
1505 # However, we need to restrict the matrix to work on the
1506 # subspace... or do we? Can't we just solve, knowing that
1507 # A(c) = u^(s+1) should have a solution in the big space,
1510 # Beware, solve_right() means that we're using COLUMN vectors.
1511 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1513 A
= u_next
.operator().matrix()
1514 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1516 # Now c is the idempotent we want, but it still lives in the subalgebra.
1517 return c
.superalgebra_element()
1522 Return my trace, the sum of my eigenvalues.
1524 In a trivial algebra, however you want to look at it, the trace is
1525 an empty sum for which we declare the result to be zero.
1529 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1536 sage: J = TrivialEJA()
1537 sage: J.zero().trace()
1541 sage: J = JordanSpinEJA(3)
1542 sage: x = sum(J.gens())
1548 sage: J = HadamardEJA(5)
1549 sage: J.one().trace()
1554 The trace of an element is a real number::
1556 sage: set_random_seed()
1557 sage: J = random_eja()
1558 sage: J.random_element().trace() in RLF
1561 The trace is linear::
1563 sage: set_random_seed()
1564 sage: J = random_eja()
1565 sage: x,y = J.random_elements(2)
1566 sage: alpha = J.base_ring().random_element()
1567 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1575 # Special case for the trivial algebra where
1576 # the trace is an empty sum.
1577 return P
.base_ring().zero()
1579 p
= P
._charpoly
_coefficients
()[r
-1]
1580 # The _charpoly_coeff function already adds the factor of
1581 # -1 to ensure that _charpoly_coeff(r-1) is really what
1582 # appears in front of t^{r-1} in the charpoly. However,
1583 # we want the negative of THAT for the trace.
1584 return -p(*self
.to_vector())
1587 def trace_inner_product(self
, other
):
1589 Return the trace inner product of myself and ``other``.
1593 sage: from mjo.eja.eja_algebra import random_eja
1597 The trace inner product is commutative, bilinear, and associative::
1599 sage: set_random_seed()
1600 sage: J = random_eja()
1601 sage: x,y,z = J.random_elements(3)
1603 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1606 sage: a = J.base_ring().random_element();
1607 sage: actual = (a*(x+z)).trace_inner_product(y)
1608 sage: expected = ( a*x.trace_inner_product(y) +
1609 ....: a*z.trace_inner_product(y) )
1610 sage: actual == expected
1612 sage: actual = x.trace_inner_product(a*(y+z))
1613 sage: expected = ( a*x.trace_inner_product(y) +
1614 ....: a*x.trace_inner_product(z) )
1615 sage: actual == expected
1618 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1622 if not other
in self
.parent():
1623 raise TypeError("'other' must live in the same algebra")
1625 return (self
*other
).trace()
1628 def trace_norm(self
):
1630 The norm of this element with respect to :meth:`trace_inner_product`.
1634 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1639 sage: J = HadamardEJA(2)
1640 sage: x = sum(J.gens())
1641 sage: x.trace_norm()
1646 sage: J = JordanSpinEJA(4)
1647 sage: x = sum(J.gens())
1648 sage: x.trace_norm()
1652 return self
.trace_inner_product(self
).sqrt()