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: x = random_eja().random_element()
50 A few examples of power-associativity::
52 sage: x = random_eja().random_element()
53 sage: x*(x*x)*(x*x) == x^5
55 sage: (x*x)*(x*x*x) == x^5
58 We also know that powers operator-commute (Koecher, Chapter
61 sage: x = random_eja().random_element()
62 sage: m = ZZ.random_element(0,10)
63 sage: n = ZZ.random_element(0,10)
64 sage: Lxm = (x^m).operator()
65 sage: Lxn = (x^n).operator()
66 sage: Lxm*Lxn == Lxn*Lxm
71 return self
.parent().one()
75 return (self
**(n
-1))*self
78 def apply_univariate_polynomial(self
, p
):
80 Apply the univariate polynomial ``p`` to this element.
82 A priori, SageMath won't allow us to apply a univariate
83 polynomial to an element of an EJA, because we don't know
84 that EJAs are rings (they are usually not associative). Of
85 course, we know that EJAs are power-associative, so the
86 operation is ultimately kosher. This function sidesteps
87 the CAS to get the answer we want and expect.
91 sage: from mjo.eja.eja_algebra import (HadamardEJA,
96 sage: R = PolynomialRing(QQ, 't')
98 sage: p = t^4 - t^3 + 5*t - 2
99 sage: J = HadamardEJA(5)
100 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
105 We should always get back an element of the algebra::
107 sage: p = PolynomialRing(AA, 't').random_element()
108 sage: J = random_eja()
109 sage: x = J.random_element()
110 sage: x.apply_univariate_polynomial(p) in J
114 if len(p
.variables()) > 1:
115 raise ValueError("not a univariate polynomial")
118 # Convert the coeficcients to the parent's base ring,
119 # because a priori they might live in an (unnecessarily)
120 # larger ring for which P.sum() would fail below.
121 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
122 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
125 def characteristic_polynomial(self
):
127 Return the characteristic polynomial of this element.
131 sage: from mjo.eja.eja_algebra import (random_eja,
136 The rank of `R^3` is three, and the minimal polynomial of
137 the identity element is `(t-1)` from which it follows that
138 the characteristic polynomial should be `(t-1)^3`::
140 sage: J = HadamardEJA(3)
141 sage: J.one().characteristic_polynomial()
142 t^3 - 3*t^2 + 3*t - 1
144 Likewise, the characteristic of the zero element in the
145 rank-three algebra `R^{n}` should be `t^{3}`::
147 sage: J = HadamardEJA(3)
148 sage: J.zero().characteristic_polynomial()
153 The characteristic polynomial of an element should evaluate
154 to zero on that element::
156 sage: x = random_eja().random_element()
157 sage: p = x.characteristic_polynomial()
158 sage: x.apply_univariate_polynomial(p).is_zero()
161 The characteristic polynomials of the zero and unit elements
162 should be what we think they are in a subalgebra, too::
164 sage: J = HadamardEJA(3)
165 sage: p1 = J.one().characteristic_polynomial()
166 sage: q1 = J.zero().characteristic_polynomial()
167 sage: b0,b1,b2 = J.gens()
168 sage: A = (b0 + 2*b1 + 3*b2).subalgebra_generated_by() # dim 3
169 sage: p2 = A.one().characteristic_polynomial()
170 sage: q2 = A.zero().characteristic_polynomial()
177 p
= self
.parent().characteristic_polynomial_of()
178 return p(*self
.to_vector())
181 def inner_product(self
, other
):
183 Return the parent algebra's inner product of myself and ``other``.
187 sage: from mjo.eja.eja_algebra import (
188 ....: ComplexHermitianEJA,
190 ....: QuaternionHermitianEJA,
191 ....: RealSymmetricEJA,
196 The inner product in the Jordan spin algebra is the usual
197 inner product on `R^n` (this example only works because the
198 basis for the Jordan algebra is the standard basis in `R^n`)::
200 sage: J = JordanSpinEJA(3)
201 sage: x = vector(QQ,[1,2,3])
202 sage: y = vector(QQ,[4,5,6])
203 sage: x.inner_product(y)
205 sage: J.from_vector(x).inner_product(J.from_vector(y))
208 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
209 multiplication is the usual matrix multiplication in `S^n`,
210 so the inner product of the identity matrix with itself
213 sage: J = RealSymmetricEJA(3)
214 sage: J.one().inner_product(J.one())
217 Likewise, the inner product on `C^n` is `<X,Y> =
218 Re(trace(X*Y))`, where we must necessarily take the real
219 part because the product of Hermitian matrices may not be
222 sage: J = ComplexHermitianEJA(3)
223 sage: J.one().inner_product(J.one())
226 Ditto for the quaternions::
228 sage: J = QuaternionHermitianEJA(2)
229 sage: J.one().inner_product(J.one())
234 Ensure that we can always compute an inner product, and that
235 it gives us back a real number::
237 sage: J = random_eja()
238 sage: x,y = J.random_elements(2)
239 sage: x.inner_product(y) in RLF
245 raise TypeError("'other' must live in the same algebra")
247 return P
.inner_product(self
, other
)
250 def operator_commutes_with(self
, other
):
252 Return whether or not this element operator-commutes
257 sage: from mjo.eja.eja_algebra import random_eja
261 The definition of a Jordan algebra says that any element
262 operator-commutes with its square::
264 sage: x = random_eja().random_element()
265 sage: x.operator_commutes_with(x^2)
270 Test Lemma 1 from Chapter III of Koecher::
272 sage: u,v = random_eja().random_elements(2)
273 sage: lhs = u.operator_commutes_with(u*v)
274 sage: rhs = v.operator_commutes_with(u^2)
278 Test the first polarization identity from my notes, Koecher
279 Chapter III, or from Baes (2.3)::
281 sage: x,y = random_eja().random_elements(2)
282 sage: Lx = x.operator()
283 sage: Ly = y.operator()
284 sage: Lxx = (x*x).operator()
285 sage: Lxy = (x*y).operator()
286 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
289 Test the second polarization identity from my notes or from
292 sage: x,y,z = random_eja().random_elements(3) # long time
293 sage: Lx = x.operator() # long time
294 sage: Ly = y.operator() # long time
295 sage: Lz = z.operator() # long time
296 sage: Lzy = (z*y).operator() # long time
297 sage: Lxy = (x*y).operator() # long time
298 sage: Lxz = (x*z).operator() # long time
299 sage: lhs = Lx*Lzy + Lz*Lxy + Ly*Lxz # long time
300 sage: rhs = Lzy*Lx + Lxy*Lz + Lxz*Ly # long time
301 sage: bool(lhs == rhs) # long time
304 Test the third polarization identity from my notes or from
307 sage: u,y,z = random_eja().random_elements(3) # long time
308 sage: Lu = u.operator() # long time
309 sage: Ly = y.operator() # long time
310 sage: Lz = z.operator() # long time
311 sage: Lzy = (z*y).operator() # long time
312 sage: Luy = (u*y).operator() # long time
313 sage: Luz = (u*z).operator() # long time
314 sage: Luyz = (u*(y*z)).operator() # long time
315 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz # long time
316 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly # long time
317 sage: bool(lhs == rhs) # long time
321 if not other
in self
.parent():
322 raise TypeError("'other' must live in the same algebra")
331 Return my determinant, the product of my eigenvalues.
335 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
337 ....: RealSymmetricEJA,
338 ....: ComplexHermitianEJA,
343 sage: J = JordanSpinEJA(2)
344 sage: x = sum( J.gens() )
350 sage: J = JordanSpinEJA(3)
351 sage: x = sum( J.gens() )
355 The determinant of the sole element in the rank-zero trivial
356 algebra is ``1``, by three paths of reasoning. First, its
357 characteristic polynomial is a constant ``1``, so the constant
358 term in that polynomial is ``1``. Second, the characteristic
359 polynomial evaluated at zero is again ``1``. And finally, the
360 (empty) product of its eigenvalues is likewise just unity::
362 sage: J = TrivialEJA()
368 An element is invertible if and only if its determinant is
371 sage: x = random_eja().random_element()
372 sage: x.is_invertible() == (x.det() != 0)
375 Ensure that the determinant is multiplicative on an associative
376 subalgebra as in Faraut and Korányi's Proposition II.2.2::
378 sage: J = random_eja().random_element().subalgebra_generated_by()
379 sage: x,y = J.random_elements(2)
380 sage: (x*y).det() == x.det()*y.det()
383 The determinant in real matrix algebras is the usual determinant::
385 sage: X = matrix.random(QQ,3)
387 sage: J1 = RealSymmetricEJA(3)
388 sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
389 sage: expected = X.det()
390 sage: actual1 = J1(X).det()
391 sage: actual2 = J2(X).det()
392 sage: actual1 == expected
394 sage: actual2 == expected
402 # Special case, since we don't get the a0=1
403 # coefficient when the rank of the algebra
405 return P
.base_ring().one()
407 p
= P
._charpoly
_coefficients
()[0]
408 # The _charpoly_coeff function already adds the factor of -1
409 # to ensure that _charpoly_coefficients()[0] is really what
410 # appears in front of t^{0} in the charpoly. However, we want
411 # (-1)^r times THAT for the determinant.
412 return ((-1)**r
)*p(*self
.to_vector())
418 Return the Jordan-multiplicative inverse of this element.
422 In general we appeal to the quadratic representation as in
423 Koecher's Theorem 12 in Chapter III, Section 5. But if the
424 parent algebra's "characteristic polynomial of" coefficients
425 happen to be cached, then we use Proposition II.2.4 in Faraut
426 and Korányi which gives a formula for the inverse based on the
427 characteristic polynomial and the Cayley-Hamilton theorem for
428 Euclidean Jordan algebras::
432 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
438 The inverse in the spin factor algebra is given in Alizadeh's
441 sage: J = JordanSpinEJA.random_instance()
442 sage: x = J.random_element()
443 sage: while not x.is_invertible():
444 ....: x = J.random_element()
445 sage: x_vec = x.to_vector()
447 sage: x_bar = x_vec[1:]
448 sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
449 sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
450 sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
451 sage: x.inverse() == J.from_vector(x_inverse)
454 Trying to invert a non-invertible element throws an error:
456 sage: JordanSpinEJA(3).zero().inverse()
457 Traceback (most recent call last):
459 ZeroDivisionError: element is not invertible
463 The identity element is its own inverse::
465 sage: J = random_eja()
466 sage: J.one().inverse() == J.one()
469 If an element has an inverse, it acts like one::
471 sage: J = random_eja()
472 sage: x = J.random_element()
473 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
476 The inverse of the inverse is what we started with::
478 sage: J = random_eja()
479 sage: x = J.random_element()
480 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
483 Proposition II.2.3 in Faraut and Korányi says that the inverse
484 of an element is the inverse of its left-multiplication operator
485 applied to the algebra's identity, when that inverse exists::
487 sage: J = random_eja()
488 sage: x = J.random_element()
489 sage: (not x.operator().is_invertible()) or (
490 ....: x.operator().inverse()(J.one()) == x.inverse() )
493 Check that the fast (cached) and slow algorithms give the same
496 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
497 sage: x = J.random_element() # long time
498 sage: while not x.is_invertible(): # long time
499 ....: x = J.random_element() # long time
500 sage: slow = x.inverse() # long time
501 sage: _ = J._charpoly_coefficients() # long time
502 sage: fast = x.inverse() # long time
503 sage: slow == fast # long time
506 not_invertible_msg
= "element is not invertible"
507 if self
.parent()._charpoly
_coefficients
.is_in_cache():
508 # We can invert using our charpoly if it will be fast to
509 # compute. If the coefficients are cached, our rank had
511 if self
.det().is_zero():
512 raise ZeroDivisionError(not_invertible_msg
)
513 r
= self
.parent().rank()
514 a
= self
.characteristic_polynomial().coefficients(sparse
=False)
515 return (-1)**(r
+1)*sum(a
[i
+1]*self
**i
for i
in range(r
))/self
.det()
518 inv
= (~self
.quadratic_representation())(self
)
519 self
.is_invertible
.set_cache(True)
521 except ZeroDivisionError:
522 self
.is_invertible
.set_cache(False)
523 raise ZeroDivisionError(not_invertible_msg
)
527 def is_invertible(self
):
529 Return whether or not this element is invertible.
533 If computing my determinant will be fast, we do so and compare
534 with zero (Proposition II.2.4 in Faraut and
535 Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
536 reduces the problem to the invertibility of my quadratic
541 sage: from mjo.eja.eja_algebra import random_eja
545 The identity element is always invertible::
547 sage: J = random_eja()
548 sage: J.one().is_invertible()
551 The zero element is never invertible in a non-trivial algebra::
553 sage: J = random_eja()
554 sage: (not J.is_trivial()) and J.zero().is_invertible()
557 Test that the fast (cached) and slow algorithms give the same
560 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
561 sage: x = J.random_element() # long time
562 sage: slow = x.is_invertible() # long time
563 sage: _ = J._charpoly_coefficients() # long time
564 sage: fast = x.is_invertible() # long time
565 sage: slow == fast # long time
569 if self
.parent().is_trivial():
574 if self
.parent()._charpoly
_coefficients
.is_in_cache():
575 # The determinant will be quicker than inverting the
576 # quadratic representation, most likely.
577 return (not self
.det().is_zero())
579 # The easiest way to determine if I'm invertible is to try.
581 inv
= (~self
.quadratic_representation())(self
)
582 self
.inverse
.set_cache(inv
)
584 except ZeroDivisionError:
588 def is_primitive_idempotent(self
):
590 Return whether or not this element is a primitive (or minimal)
593 A primitive idempotent is a non-zero idempotent that is not
594 the sum of two other non-zero idempotents. Remark 2.7.15 in
595 Baes shows that this is what he refers to as a "minimal
598 An element of a Euclidean Jordan algebra is a minimal idempotent
599 if it :meth:`is_idempotent` and if its Peirce subalgebra
600 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
605 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
606 ....: RealSymmetricEJA,
612 This method is sloooooow.
616 The spectral decomposition of a non-regular element should always
617 contain at least one non-minimal idempotent::
619 sage: J = RealSymmetricEJA(3)
620 sage: x = sum(J.gens())
623 sage: [ c.is_primitive_idempotent()
624 ....: for (l,c) in x.spectral_decomposition() ]
627 On the other hand, the spectral decomposition of a regular
628 element should always be in terms of minimal idempotents::
630 sage: J = JordanSpinEJA(4)
631 sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
634 sage: [ c.is_primitive_idempotent()
635 ....: for (l,c) in x.spectral_decomposition() ]
640 The identity element is minimal only in an EJA of rank one::
642 sage: J = random_eja()
643 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
646 A non-idempotent cannot be a minimal idempotent::
648 sage: J = JordanSpinEJA(4)
649 sage: x = J.random_element()
650 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
653 Proposition 2.7.19 in Baes says that an element is a minimal
654 idempotent if and only if it's idempotent with trace equal to
657 sage: J = JordanSpinEJA(4)
658 sage: x = J.random_element()
659 sage: expected = (x.is_idempotent() and x.trace() == 1)
660 sage: actual = x.is_primitive_idempotent()
661 sage: actual == expected
664 Primitive idempotents must be non-zero::
666 sage: J = random_eja()
667 sage: J.zero().is_idempotent()
669 sage: J.zero().is_primitive_idempotent()
672 As a consequence of the fact that primitive idempotents must
673 be non-zero, there are no primitive idempotents in a trivial
674 Euclidean Jordan algebra::
676 sage: J = TrivialEJA()
677 sage: J.one().is_idempotent()
679 sage: J.one().is_primitive_idempotent()
683 if not self
.is_idempotent():
689 (_
,_
,J1
) = self
.parent().peirce_decomposition(self
)
690 return (J1
.dimension() == 1)
693 def is_nilpotent(self
):
695 Return whether or not some power of this element is zero.
699 We use Theorem 5 in Chapter III of Koecher, which says that
700 an element ``x`` is nilpotent if and only if ``x.operator()``
701 is nilpotent. And it is a basic fact of linear algebra that
702 an operator on an `n`-dimensional space is nilpotent if and
703 only if, when raised to the `n`th power, it equals the zero
704 operator (for example, see Axler Corollary 8.8).
708 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
713 sage: J = JordanSpinEJA(3)
714 sage: x = sum(J.gens())
715 sage: x.is_nilpotent()
720 The identity element is never nilpotent, except in a trivial EJA::
722 sage: J = random_eja()
723 sage: J.one().is_nilpotent() and not J.is_trivial()
726 The additive identity is always nilpotent::
728 sage: random_eja().zero().is_nilpotent()
733 zero_operator
= P
.zero().operator()
734 return self
.operator()**P
.dimension() == zero_operator
737 def is_regular(self
):
739 Return whether or not this is a regular element.
743 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
748 The identity element always has degree one, but any element
749 linearly-independent from it is regular::
751 sage: J = JordanSpinEJA(5)
752 sage: J.one().is_regular()
754 sage: b0, b1, b2, b3, b4 = J.gens()
757 sage: for x in J.gens():
758 ....: (J.one() + x).is_regular()
767 The zero element should never be regular, unless the parent
768 algebra has dimension less than or equal to one::
770 sage: J = random_eja()
771 sage: J.dimension() <= 1 or not J.zero().is_regular()
774 The unit element isn't regular unless the algebra happens to
775 consist of only its scalar multiples::
777 sage: J = random_eja()
778 sage: J.dimension() <= 1 or not J.one().is_regular()
782 return self
.degree() == self
.parent().rank()
787 Return the degree of this element, which is defined to be
788 the degree of its minimal polynomial.
796 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
801 sage: J = JordanSpinEJA(4)
802 sage: J.one().degree()
804 sage: b0,b1,b2,b3 = J.gens()
805 sage: (b0 - b1).degree()
808 In the spin factor algebra (of rank two), all elements that
809 aren't multiples of the identity are regular::
811 sage: J = JordanSpinEJA.random_instance()
812 sage: n = J.dimension()
813 sage: x = J.random_element()
814 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
819 The zero and unit elements are both of degree one in nontrivial
822 sage: J = random_eja()
823 sage: d = J.zero().degree()
824 sage: (J.is_trivial() and d == 0) or d == 1
826 sage: d = J.one().degree()
827 sage: (J.is_trivial() and d == 0) or d == 1
830 Our implementation agrees with the definition::
832 sage: x = random_eja().random_element()
833 sage: x.degree() == x.minimal_polynomial().degree()
837 n
= self
.parent().dimension()
840 # The minimal polynomial is an empty product, i.e. the
841 # constant polynomial "1" having degree zero.
844 # The minimal polynomial of zero in a nontrivial algebra
845 # is "t", and is of degree one.
848 # If this is a nonzero element of a nontrivial algebra, it
849 # has degree at least one. It follows that, in an algebra
850 # of dimension one, the degree must be actually one.
853 # BEWARE: The subalgebra_generated_by() method uses the result
854 # of this method to construct a basis for the subalgebra. That
855 # means, in particular, that we cannot implement this method
856 # as ``self.subalgebra_generated_by().dimension()``.
858 # Algorithm: keep appending (vector representations of) powers
859 # self as rows to a matrix and echelonizing it. When its rank
860 # stops increasing, we've reached a redundancy.
862 # Given the special cases above, we can assume that "self" is
863 # nonzero, the algebra is nontrivial, and that its dimension
865 M
= matrix([(self
.parent().one()).to_vector()])
868 # Specifying the row-reduction algorithm can e.g. help over
869 # AA because it avoids the RecursionError that gets thrown
870 # when we have to look too hard for a root.
872 # Beware: QQ supports an entirely different set of "algorithm"
873 # keywords than do AA and RR.
875 from sage
.rings
.all
import QQ
876 if self
.parent().base_ring() is not QQ
:
877 algo
= "scaled_partial_pivoting"
880 M
= matrix(M
.rows() + [(self
**d
).to_vector()])
883 if new_rank
== old_rank
:
892 def left_matrix(self
):
894 Our parent class defines ``left_matrix`` and ``matrix``
895 methods whose names are misleading. We don't want them.
897 raise NotImplementedError("use operator().matrix() instead")
902 def minimal_polynomial(self
):
904 Return the minimal polynomial of this element,
905 as a function of the variable `t`.
909 We restrict ourselves to the associative subalgebra
910 generated by this element, and then return the minimal
911 polynomial of this element's operator matrix (in that
912 subalgebra). This works by Baes Proposition 2.3.16.
916 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
917 ....: RealSymmetricEJA,
923 Keeping in mind that the polynomial ``1`` evaluates the identity
924 element (also the zero element) of the trivial algebra, it is clear
925 that the polynomial ``1`` is the minimal polynomial of the only
926 element in a trivial algebra::
928 sage: J = TrivialEJA()
929 sage: J.one().minimal_polynomial()
931 sage: J.zero().minimal_polynomial()
936 The minimal polynomial of the identity and zero elements are
937 always the same, except in trivial algebras where the minimal
938 polynomial of the unit/zero element is ``1``::
940 sage: J = random_eja()
941 sage: mu = J.one().minimal_polynomial()
942 sage: t = mu.parent().gen()
943 sage: mu + int(J.is_trivial())*(t-2)
945 sage: mu = J.zero().minimal_polynomial()
946 sage: t = mu.parent().gen()
947 sage: mu + int(J.is_trivial())*(t-1)
950 The degree of an element is (by one definition) the degree
951 of its minimal polynomial::
953 sage: x = random_eja().random_element()
954 sage: x.degree() == x.minimal_polynomial().degree()
957 The minimal polynomial and the characteristic polynomial coincide
958 and are known (see Alizadeh, Example 11.11) for all elements of
959 the spin factor algebra that aren't scalar multiples of the
960 identity. We require the dimension of the algebra to be at least
961 two here so that said elements actually exist::
963 sage: d_max = JordanSpinEJA._max_random_instance_dimension()
964 sage: n = ZZ.random_element(2, max(2,d_max))
965 sage: J = JordanSpinEJA(n)
966 sage: y = J.random_element()
967 sage: while y == y.coefficient(0)*J.one():
968 ....: y = J.random_element()
969 sage: y0 = y.to_vector()[0]
970 sage: y_bar = y.to_vector()[1:]
971 sage: actual = y.minimal_polynomial()
972 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
973 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
974 sage: bool(actual == expected)
977 The minimal polynomial should always kill its element::
979 sage: x = random_eja().random_element()
980 sage: p = x.minimal_polynomial()
981 sage: x.apply_univariate_polynomial(p)
984 The minimal polynomial is invariant under a change of basis,
985 and in particular, a re-scaling of the basis::
987 sage: d_max = RealSymmetricEJA._max_random_instance_dimension()
988 sage: d = ZZ.random_element(1, d_max)
989 sage: n = RealSymmetricEJA._max_random_instance_size(d)
990 sage: J1 = RealSymmetricEJA(n)
991 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
992 sage: X = random_matrix(AA,n)
993 sage: X = X*X.transpose()
996 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
1001 # Pretty sure we know what the minimal polynomial of
1002 # the zero operator is going to be. This ensures
1003 # consistency of e.g. the polynomial variable returned
1004 # in the "normal" case without us having to think about it.
1005 return self
.operator().minimal_polynomial()
1007 # If we don't orthonormalize the subalgebra's basis, then the
1008 # first two monomials in the subalgebra will be self^0 and
1009 # self^1... assuming that self^1 is not a scalar multiple of
1010 # self^0 (the unit element). We special case these to avoid
1011 # having to solve a system to coerce self into the subalgebra.
1012 A
= self
.subalgebra_generated_by(orthonormalize
=False)
1014 if A
.dimension() == 1:
1015 # Does a solve to find the scalar multiple alpha such that
1016 # alpha*unit = self. We have to do this because the basis
1017 # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
1018 unit
= self
.parent().one()
1019 alpha
= self
.to_vector() / unit
.to_vector()
1020 return (unit
.operator()*alpha
).minimal_polynomial()
1022 # If the dimension of the subalgebra is >= 2, then we just
1023 # use the second basis element.
1024 return A
.monomial(1).operator().minimal_polynomial()
1028 def to_matrix(self
):
1030 Return an (often more natural) representation of this element as a
1033 Every finite-dimensional Euclidean Jordan Algebra is a direct
1034 sum of five simple algebras, four of which comprise Hermitian
1035 matrices. This method returns a "natural" matrix
1036 representation of this element as either a Hermitian matrix or
1041 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1043 ....: QuaternionHermitianEJA,
1044 ....: RealSymmetricEJA)
1048 sage: J = ComplexHermitianEJA(3)
1051 sage: J.one().to_matrix()
1062 sage: J = QuaternionHermitianEJA(2)
1065 sage: J.one().to_matrix()
1072 This also works in Cartesian product algebras::
1074 sage: J1 = HadamardEJA(1)
1075 sage: J2 = RealSymmetricEJA(2)
1076 sage: J = cartesian_product([J1,J2])
1077 sage: x = sum(J.gens())
1078 sage: x.to_matrix()[0]
1080 sage: x.to_matrix()[1]
1081 [ 1 0.7071067811865475?]
1082 [0.7071067811865475? 1]
1085 B
= self
.parent().matrix_basis()
1086 W
= self
.parent().matrix_space()
1088 # This is just a manual "from_vector()", but of course
1089 # matrix spaces aren't vector spaces in sage, so they
1090 # don't have a from_vector() method.
1091 return W
.linear_combination( zip(B
, self
.to_vector()) )
1097 The norm of this element with respect to :meth:`inner_product`.
1101 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1106 sage: J = HadamardEJA(2)
1107 sage: x = sum(J.gens())
1113 sage: J = JordanSpinEJA(4)
1114 sage: x = sum(J.gens())
1119 return self
.inner_product(self
).sqrt()
1124 Return the left-multiplication-by-this-element
1125 operator on the ambient algebra.
1129 sage: from mjo.eja.eja_algebra import random_eja
1133 sage: J = random_eja()
1134 sage: x,y = J.random_elements(2)
1135 sage: x.operator()(y) == x*y
1137 sage: y.operator()(x) == x*y
1142 left_mult_by_self
= lambda y
: self
*y
1143 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1144 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1147 def quadratic_representation(self
, other
=None):
1149 Return the quadratic representation of this element.
1153 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1158 The explicit form in the spin factor algebra is given by
1159 Alizadeh's Example 11.12::
1161 sage: x = JordanSpinEJA.random_instance().random_element()
1162 sage: x_vec = x.to_vector()
1163 sage: Q = matrix.identity(x.base_ring(), 0)
1164 sage: n = x_vec.degree()
1167 ....: x_bar = x_vec[1:]
1168 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1169 ....: B = 2*x0*x_bar.row()
1170 ....: C = 2*x0*x_bar.column()
1171 ....: D = matrix.identity(x.base_ring(), n-1)
1172 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1173 ....: D = D + 2*x_bar.tensor_product(x_bar)
1174 ....: Q = matrix.block(2,2,[A,B,C,D])
1175 sage: Q == x.quadratic_representation().matrix()
1178 Test all of the properties from Theorem 11.2 in Alizadeh::
1180 sage: J = random_eja()
1181 sage: x,y = J.random_elements(2)
1182 sage: Lx = x.operator()
1183 sage: Lxx = (x*x).operator()
1184 sage: Qx = x.quadratic_representation()
1185 sage: Qy = y.quadratic_representation()
1186 sage: Qxy = x.quadratic_representation(y)
1187 sage: Qex = J.one().quadratic_representation(x)
1188 sage: n = ZZ.random_element(10)
1189 sage: Qxn = (x^n).quadratic_representation()
1193 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1196 Property 2 (multiply on the right for :trac:`28272`):
1198 sage: alpha = J.base_ring().random_element()
1199 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1204 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1207 sage: not x.is_invertible() or (
1210 ....: x.inverse().quadratic_representation() )
1213 sage: Qxy(J.one()) == x*y
1218 sage: not x.is_invertible() or (
1219 ....: x.quadratic_representation(x.inverse())*Qx
1220 ....: == Qx*x.quadratic_representation(x.inverse()) )
1223 sage: not x.is_invertible() or (
1224 ....: x.quadratic_representation(x.inverse())*Qx
1226 ....: 2*Lx*Qex - Qx )
1229 sage: 2*Lx*Qex - Qx == Lxx
1234 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1244 sage: not x.is_invertible() or (
1245 ....: Qx*x.inverse().operator() == Lx )
1250 sage: not x.operator_commutes_with(y) or (
1251 ....: Qx(y)^n == Qxn(y^n) )
1257 elif not other
in self
.parent():
1258 raise TypeError("'other' must live in the same algebra")
1261 M
= other
.operator()
1262 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1266 def spectral_decomposition(self
):
1268 Return the unique spectral decomposition of this element.
1272 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1273 element's left-multiplication-by operator to the subalgebra it
1274 generates. We then compute the spectral decomposition of that
1275 operator, and the spectral projectors we get back must be the
1276 left-multiplication-by operators for the idempotents we
1277 seek. Thus applying them to the identity element gives us those
1280 Since the eigenvalues are required to be distinct, we take
1281 the spectral decomposition of the zero element to be zero
1282 times the identity element of the algebra (which is idempotent,
1287 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1291 The spectral decomposition of the identity is ``1`` times itself,
1292 and the spectral decomposition of zero is ``0`` times the identity::
1294 sage: J = RealSymmetricEJA(3)
1297 sage: J.one().spectral_decomposition()
1299 sage: J.zero().spectral_decomposition()
1304 sage: J = RealSymmetricEJA(4)
1305 sage: x = sum(J.gens())
1306 sage: sd = x.spectral_decomposition()
1311 sage: c0.inner_product(c1) == 0
1313 sage: c0.is_idempotent()
1315 sage: c1.is_idempotent()
1317 sage: c0 + c1 == J.one()
1319 sage: l0*c0 + l1*c1 == x
1322 The spectral decomposition should work in subalgebras, too::
1324 sage: J = RealSymmetricEJA(4)
1325 sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
1326 sage: A = 2*b5 - 2*b8
1327 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1328 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1329 sage: (f0, f1, f2) = J1.gens()
1330 sage: f0.spectral_decomposition()
1331 [(0, 1.000000000000000?*c2), (1, 1.000000000000000?*c0)]
1334 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1336 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1337 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1340 def subalgebra_generated_by(self
, **kwargs
):
1342 Return the associative subalgebra of the parent EJA generated
1345 Since our parent algebra is unital, we want "subalgebra" to mean
1346 "unital subalgebra" as well; thus the subalgebra that an element
1347 generates will itself be a Euclidean Jordan algebra after
1348 restricting the algebra operations appropriately. This is the
1349 subalgebra that Faraut and Korányi work with in section II.2, for
1354 sage: from mjo.eja.eja_algebra import (random_eja,
1356 ....: RealSymmetricEJA)
1360 We can create subalgebras of Cartesian product EJAs that are not
1361 themselves Cartesian product EJAs (they're just "regular" EJAs)::
1363 sage: J1 = HadamardEJA(3)
1364 sage: J2 = RealSymmetricEJA(2)
1365 sage: J = cartesian_product([J1,J2])
1366 sage: J.one().subalgebra_generated_by()
1367 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
1371 This subalgebra, being composed of only powers, is associative::
1373 sage: x0 = random_eja().random_element()
1374 sage: A = x0.subalgebra_generated_by()
1375 sage: x,y,z = A.random_elements(3)
1376 sage: (x*y)*z == x*(y*z)
1379 Squaring in the subalgebra should work the same as in
1382 sage: x = random_eja().random_element()
1383 sage: A = x.subalgebra_generated_by()
1384 sage: A(x^2) == A(x)*A(x)
1387 By definition, the subalgebra generated by the zero element is
1388 the one-dimensional algebra generated by the identity
1389 element... unless the original algebra was trivial, in which
1390 case the subalgebra is trivial too::
1392 sage: A = random_eja().zero().subalgebra_generated_by()
1393 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1397 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1398 A
= self
.parent().subalgebra(powers
,
1403 A
.one
.set_cache(A(self
.parent().one()))
1407 def subalgebra_idempotent(self
):
1409 Find an idempotent in the associative subalgebra I generate
1410 using Proposition 2.3.5 in Baes.
1414 sage: from mjo.eja.eja_algebra import random_eja
1418 Ensure that we can find an idempotent in a non-trivial algebra
1419 where there are non-nilpotent elements, or that we get the dumb
1420 solution in the trivial algebra::
1422 sage: J = random_eja()
1423 sage: x = J.random_element()
1424 sage: while x.is_nilpotent() and not J.is_trivial():
1425 ....: x = J.random_element()
1426 sage: c = x.subalgebra_idempotent()
1431 if self
.parent().is_trivial():
1434 if self
.is_nilpotent():
1435 raise ValueError("this only works with non-nilpotent elements!")
1437 J
= self
.subalgebra_generated_by()
1440 # The image of the matrix of left-u^m-multiplication
1441 # will be minimal for some natural number s...
1443 minimal_dim
= J
.dimension()
1444 for i
in range(1, minimal_dim
):
1445 this_dim
= (u
**i
).operator().matrix().image().dimension()
1446 if this_dim
< minimal_dim
:
1447 minimal_dim
= this_dim
1450 # Now minimal_matrix should correspond to the smallest
1451 # non-zero subspace in Baes's (or really, Koecher's)
1454 # However, we need to restrict the matrix to work on the
1455 # subspace... or do we? Can't we just solve, knowing that
1456 # A(c) = u^(s+1) should have a solution in the big space,
1459 # Beware, solve_right() means that we're using COLUMN vectors.
1460 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1462 A
= u_next
.operator().matrix()
1463 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1465 # Now c is the idempotent we want, but it still lives in the subalgebra.
1466 return c
.superalgebra_element()
1471 Return my trace, the sum of my eigenvalues.
1473 In a trivial algebra, however you want to look at it, the trace is
1474 an empty sum for which we declare the result to be zero.
1478 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1485 sage: J = TrivialEJA()
1486 sage: J.zero().trace()
1490 sage: J = JordanSpinEJA(3)
1491 sage: x = sum(J.gens())
1497 sage: J = HadamardEJA(5)
1498 sage: J.one().trace()
1503 The trace of an element is a real number::
1505 sage: J = random_eja()
1506 sage: J.random_element().trace() in RLF
1509 The trace is linear::
1511 sage: J = random_eja()
1512 sage: x,y = J.random_elements(2)
1513 sage: alpha = J.base_ring().random_element()
1514 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1517 The trace of a square is nonnegative::
1519 sage: x = random_eja().random_element()
1520 sage: (x*x).trace() >= 0
1528 # Special case for the trivial algebra where
1529 # the trace is an empty sum.
1530 return P
.base_ring().zero()
1532 p
= P
._charpoly
_coefficients
()[r
-1]
1533 # The _charpoly_coeff function already adds the factor of
1534 # -1 to ensure that _charpoly_coeff(r-1) is really what
1535 # appears in front of t^{r-1} in the charpoly. However,
1536 # we want the negative of THAT for the trace.
1537 return -p(*self
.to_vector())
1540 def trace_inner_product(self
, other
):
1542 Return the trace inner product of myself and ``other``.
1546 sage: from mjo.eja.eja_algebra import random_eja
1550 The trace inner product is commutative, bilinear, and associative::
1552 sage: J = random_eja()
1553 sage: x,y,z = J.random_elements(3)
1555 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1558 sage: a = J.base_ring().random_element();
1559 sage: actual = (a*(x+z)).trace_inner_product(y)
1560 sage: expected = ( a*x.trace_inner_product(y) +
1561 ....: a*z.trace_inner_product(y) )
1562 sage: actual == expected
1564 sage: actual = x.trace_inner_product(a*(y+z))
1565 sage: expected = ( a*x.trace_inner_product(y) +
1566 ....: a*x.trace_inner_product(z) )
1567 sage: actual == expected
1570 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1574 if not other
in self
.parent():
1575 raise TypeError("'other' must live in the same algebra")
1577 return (self
*other
).trace()
1580 def trace_norm(self
):
1582 The norm of this element with respect to :meth:`trace_inner_product`.
1586 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1591 sage: J = HadamardEJA(2)
1592 sage: x = sum(J.gens())
1593 sage: x.trace_norm()
1598 sage: J = JordanSpinEJA(4)
1599 sage: x = sum(J.gens())
1600 sage: x.trace_norm()
1604 return self
.trace_inner_product(self
).sqrt()
1607 class CartesianProductEJAElement(FiniteDimensionalEJAElement
):
1610 Compute the determinant of this product-element using the
1611 determianants of its factors.
1613 This result Follows from the spectral decomposition of (say)
1614 the pair `(x,y)` in terms of the Jordan frame `\left\{ (c_1,
1615 0),(c_2, 0),...,(0,d_1),(0,d_2),... \right\}.
1617 from sage
.misc
.misc_c
import prod
1618 return prod( f
.det() for f
in self
.cartesian_factors() )
1620 def to_matrix(self
):
1621 # An override is necessary to call our custom _scale().
1622 B
= self
.parent().matrix_basis()
1623 W
= self
.parent().matrix_space()
1625 # Aaaaand linear combinations don't work in Cartesian
1626 # product spaces, even though they provide a method with
1627 # that name. This is hidden behind an "if" because the
1628 # _scale() function is slow.
1629 pairs
= zip(B
, self
.to_vector())
1630 return W
.sum( _scale(b
, alpha
) for (b
,alpha
) in pairs
)