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: x0 = random_eja().random_element()
379 sage: J = x0.subalgebra_generated_by(orthonormalize=False)
380 sage: x,y = J.random_elements(2)
381 sage: (x*y).det() == x.det()*y.det()
384 The determinant in real matrix algebras is the usual determinant::
386 sage: X = matrix.random(QQ,3)
388 sage: J1 = RealSymmetricEJA(3)
389 sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
390 sage: expected = X.det()
391 sage: actual1 = J1(X).det()
392 sage: actual2 = J2(X).det()
393 sage: actual1 == expected
395 sage: actual2 == expected
403 # Special case, since we don't get the a0=1
404 # coefficient when the rank of the algebra
406 return P
.base_ring().one()
408 p
= P
._charpoly
_coefficients
()[0]
409 # The _charpoly_coeff function already adds the factor of -1
410 # to ensure that _charpoly_coefficients()[0] is really what
411 # appears in front of t^{0} in the charpoly. However, we want
412 # (-1)^r times THAT for the determinant.
413 return ((-1)**r
)*p(*self
.to_vector())
419 Return the Jordan-multiplicative inverse of this element.
423 In general we appeal to the quadratic representation as in
424 Koecher's Theorem 12 in Chapter III, Section 5. But if the
425 parent algebra's "characteristic polynomial of" coefficients
426 happen to be cached, then we use Proposition II.2.4 in Faraut
427 and Korányi which gives a formula for the inverse based on the
428 characteristic polynomial and the Cayley-Hamilton theorem for
429 Euclidean Jordan algebras::
433 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
439 The inverse in the spin factor algebra is given in Alizadeh's
442 sage: J = JordanSpinEJA.random_instance()
443 sage: x = J.random_element()
444 sage: while not x.is_invertible():
445 ....: x = J.random_element()
446 sage: x_vec = x.to_vector()
448 sage: x_bar = x_vec[1:]
449 sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
450 sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
451 sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
452 sage: x.inverse() == J.from_vector(x_inverse)
455 Trying to invert a non-invertible element throws an error:
457 sage: JordanSpinEJA(3).zero().inverse()
458 Traceback (most recent call last):
460 ZeroDivisionError: element is not invertible
464 The identity element is its own inverse::
466 sage: J = random_eja()
467 sage: J.one().inverse() == J.one()
470 If an element has an inverse, it acts like one::
472 sage: J = random_eja()
473 sage: x = J.random_element()
474 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
477 The inverse of the inverse is what we started with::
479 sage: J = random_eja()
480 sage: x = J.random_element()
481 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
484 Proposition II.2.3 in Faraut and Korányi says that the inverse
485 of an element is the inverse of its left-multiplication operator
486 applied to the algebra's identity, when that inverse exists::
488 sage: J = random_eja() # long time
489 sage: x = J.random_element() # long time
490 sage: (not x.operator().is_invertible()) or ( # long time
491 ....: x.operator().inverse()(J.one()) # long time
493 ....: x.inverse() ) # long time
496 Check that the fast (cached) and slow algorithms give the same
499 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
500 sage: x = J.random_element() # long time
501 sage: while not x.is_invertible(): # long time
502 ....: x = J.random_element() # long time
503 sage: slow = x.inverse() # long time
504 sage: _ = J._charpoly_coefficients() # long time
505 sage: fast = x.inverse() # long time
506 sage: slow == fast # long time
509 not_invertible_msg
= "element is not invertible"
511 algebra
= self
.parent()
512 if algebra
._charpoly
_coefficients
.is_in_cache():
513 # We can invert using our charpoly if it will be fast to
514 # compute. If the coefficients are cached, our rank had
516 if self
.det().is_zero():
517 raise ZeroDivisionError(not_invertible_msg
)
519 a
= self
.characteristic_polynomial().coefficients(sparse
=False)
520 return (-1)**(r
+1)*algebra
.sum(a
[i
+1]*self
**i
521 for i
in range(r
))/self
.det()
524 inv
= (~self
.quadratic_representation())(self
)
525 self
.is_invertible
.set_cache(True)
527 except ZeroDivisionError:
528 self
.is_invertible
.set_cache(False)
529 raise ZeroDivisionError(not_invertible_msg
)
533 def is_invertible(self
):
535 Return whether or not this element is invertible.
539 If computing my determinant will be fast, we do so and compare
540 with zero (Proposition II.2.4 in Faraut and
541 Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
542 reduces the problem to the invertibility of my quadratic
547 sage: from mjo.eja.eja_algebra import random_eja
551 The identity element is always invertible::
553 sage: J = random_eja()
554 sage: J.one().is_invertible()
557 The zero element is never invertible in a non-trivial algebra::
559 sage: J = random_eja()
560 sage: (not J.is_trivial()) and J.zero().is_invertible()
563 Test that the fast (cached) and slow algorithms give the same
566 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
567 sage: x = J.random_element() # long time
568 sage: slow = x.is_invertible() # long time
569 sage: _ = J._charpoly_coefficients() # long time
570 sage: fast = x.is_invertible() # long time
571 sage: slow == fast # long time
575 if self
.parent().is_trivial():
580 if self
.parent()._charpoly
_coefficients
.is_in_cache():
581 # The determinant will be quicker than inverting the
582 # quadratic representation, most likely.
583 return (not self
.det().is_zero())
585 # The easiest way to determine if I'm invertible is to try.
587 inv
= (~self
.quadratic_representation())(self
)
588 self
.inverse
.set_cache(inv
)
590 except ZeroDivisionError:
594 def is_primitive_idempotent(self
):
596 Return whether or not this element is a primitive (or minimal)
599 A primitive idempotent is a non-zero idempotent that is not
600 the sum of two other non-zero idempotents. Remark 2.7.15 in
601 Baes shows that this is what he refers to as a "minimal
604 An element of a Euclidean Jordan algebra is a minimal idempotent
605 if it :meth:`is_idempotent` and if its Peirce subalgebra
606 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
611 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
612 ....: RealSymmetricEJA,
618 This method is sloooooow.
622 The spectral decomposition of a non-regular element should always
623 contain at least one non-minimal idempotent::
625 sage: J = RealSymmetricEJA(3)
626 sage: x = sum(J.gens())
629 sage: [ c.is_primitive_idempotent()
630 ....: for (l,c) in x.spectral_decomposition() ]
633 On the other hand, the spectral decomposition of a regular
634 element should always be in terms of minimal idempotents::
636 sage: J = JordanSpinEJA(4)
637 sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
640 sage: [ c.is_primitive_idempotent()
641 ....: for (l,c) in x.spectral_decomposition() ]
646 The identity element is minimal only in an EJA of rank one::
648 sage: J = random_eja()
649 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
652 A non-idempotent cannot be a minimal idempotent::
654 sage: J = JordanSpinEJA(4)
655 sage: x = J.random_element()
656 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
659 Proposition 2.7.19 in Baes says that an element is a minimal
660 idempotent if and only if it's idempotent with trace equal to
663 sage: J = JordanSpinEJA(4)
664 sage: x = J.random_element()
665 sage: expected = (x.is_idempotent() and x.trace() == 1)
666 sage: actual = x.is_primitive_idempotent()
667 sage: actual == expected
670 Primitive idempotents must be non-zero::
672 sage: J = random_eja()
673 sage: J.zero().is_idempotent()
675 sage: J.zero().is_primitive_idempotent()
678 As a consequence of the fact that primitive idempotents must
679 be non-zero, there are no primitive idempotents in a trivial
680 Euclidean Jordan algebra::
682 sage: J = TrivialEJA()
683 sage: J.one().is_idempotent()
685 sage: J.one().is_primitive_idempotent()
689 if not self
.is_idempotent():
695 (_
,_
,J1
) = self
.parent().peirce_decomposition(self
)
696 return (J1
.dimension() == 1)
699 def is_nilpotent(self
):
701 Return whether or not some power of this element is zero.
705 We use Theorem 5 in Chapter III of Koecher, which says that
706 an element ``x`` is nilpotent if and only if ``x.operator()``
707 is nilpotent. And it is a basic fact of linear algebra that
708 an operator on an `n`-dimensional space is nilpotent if and
709 only if, when raised to the `n`th power, it equals the zero
710 operator (for example, see Axler Corollary 8.8).
714 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
719 sage: J = JordanSpinEJA(3)
720 sage: x = sum(J.gens())
721 sage: x.is_nilpotent()
726 The identity element is never nilpotent, except in a trivial EJA::
728 sage: J = random_eja()
729 sage: J.one().is_nilpotent() and not J.is_trivial()
732 The additive identity is always nilpotent::
734 sage: random_eja().zero().is_nilpotent()
739 zero_operator
= P
.zero().operator()
740 return self
.operator()**P
.dimension() == zero_operator
743 def is_regular(self
):
745 Return whether or not this is a regular element.
749 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
754 The identity element always has degree one, but any element
755 linearly-independent from it is regular::
757 sage: J = JordanSpinEJA(5)
758 sage: J.one().is_regular()
760 sage: b0, b1, b2, b3, b4 = J.gens()
763 sage: for x in J.gens():
764 ....: (J.one() + x).is_regular()
773 The zero element should never be regular, unless the parent
774 algebra has dimension less than or equal to one::
776 sage: J = random_eja()
777 sage: J.dimension() <= 1 or not J.zero().is_regular()
780 The unit element isn't regular unless the algebra happens to
781 consist of only its scalar multiples::
783 sage: J = random_eja()
784 sage: J.dimension() <= 1 or not J.one().is_regular()
788 return self
.degree() == self
.parent().rank()
793 Return the degree of this element, which is defined to be
794 the degree of its minimal polynomial.
798 First we handle the special cases where the algebra is
799 trivial, this element is zero, or the dimension of the algebra
800 is one and this element is not zero. With those out of the
801 way, we may assume that ``self`` is nonzero, the algebra is
802 nontrivial, and that the dimension of the algebra is at least
805 Beginning with the algebra's unit element (power zero), we add
806 successive (basis representations of) powers of this element
807 to a matrix, row-reducing at each step. After row-reducing, we
808 check the rank of the matrix. If adding a row and row-reducing
809 does not increase the rank of the matrix at any point, the row
810 we've just added lives in the span of the previous ones; thus
811 the corresponding power of ``self`` lives in the span of its
812 lesser powers. When that happens, the degree of the minimal
813 polynomial is the rank of the matrix; if it never happens, the
814 degree must be the dimension of the entire space.
818 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
823 sage: J = JordanSpinEJA(4)
824 sage: J.one().degree()
826 sage: b0,b1,b2,b3 = J.gens()
827 sage: (b0 - b1).degree()
830 In the spin factor algebra (of rank two), all elements that
831 aren't multiples of the identity are regular::
833 sage: J = JordanSpinEJA.random_instance()
834 sage: n = J.dimension()
835 sage: x = J.random_element()
836 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
841 The zero and unit elements are both of degree one in nontrivial
844 sage: J = random_eja()
845 sage: d = J.zero().degree()
846 sage: (J.is_trivial() and d == 0) or d == 1
848 sage: d = J.one().degree()
849 sage: (J.is_trivial() and d == 0) or d == 1
852 Our implementation agrees with the definition::
854 sage: x = random_eja().random_element()
855 sage: x.degree() == x.minimal_polynomial().degree()
858 n
= self
.parent().dimension()
861 # The minimal polynomial is an empty product, i.e. the
862 # constant polynomial "1" having degree zero.
865 # The minimal polynomial of zero in a nontrivial algebra
866 # is "t", and is of degree one.
869 # If this is a nonzero element of a nontrivial algebra, it
870 # has degree at least one. It follows that, in an algebra
871 # of dimension one, the degree must be actually one.
874 # BEWARE: The subalgebra_generated_by() method uses the result
875 # of this method to construct a basis for the subalgebra. That
876 # means, in particular, that we cannot implement this method
877 # as ``self.subalgebra_generated_by().dimension()``.
879 # Algorithm: keep appending (vector representations of) powers
880 # self as rows to a matrix and echelonizing it. When its rank
881 # stops increasing, we've reached a redundancy.
883 # Given the special cases above, we can assume that "self" is
884 # nonzero, the algebra is nontrivial, and that its dimension
886 M
= matrix([(self
.parent().one()).to_vector()])
889 # Specifying the row-reduction algorithm can e.g. help over
890 # AA because it avoids the RecursionError that gets thrown
891 # when we have to look too hard for a root.
893 # Beware: QQ supports an entirely different set of "algorithm"
894 # keywords than do AA and RR.
896 from sage
.rings
.all
import QQ
897 if self
.parent().base_ring() is not QQ
:
898 algo
= "scaled_partial_pivoting"
901 M
= matrix(M
.rows() + [(self
**d
).to_vector()])
904 if new_rank
== old_rank
:
913 def left_matrix(self
):
915 Our parent class defines ``left_matrix`` and ``matrix``
916 methods whose names are misleading. We don't want them.
918 raise NotImplementedError("use operator().matrix() instead")
923 def minimal_polynomial(self
):
925 Return the minimal polynomial of this element,
926 as a function of the variable `t`.
930 We restrict ourselves to the associative subalgebra
931 generated by this element, and then return the minimal
932 polynomial of this element's operator matrix (in that
933 subalgebra). This works by Baes Proposition 2.3.16.
937 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
938 ....: RealSymmetricEJA,
944 Keeping in mind that the polynomial ``1`` evaluates the identity
945 element (also the zero element) of the trivial algebra, it is clear
946 that the polynomial ``1`` is the minimal polynomial of the only
947 element in a trivial algebra::
949 sage: J = TrivialEJA()
950 sage: J.one().minimal_polynomial()
952 sage: J.zero().minimal_polynomial()
957 The minimal polynomial of the identity and zero elements are
958 always the same, except in trivial algebras where the minimal
959 polynomial of the unit/zero element is ``1``::
961 sage: J = random_eja()
962 sage: mu = J.one().minimal_polynomial()
963 sage: t = mu.parent().gen()
964 sage: mu + int(J.is_trivial())*(t-2)
966 sage: mu = J.zero().minimal_polynomial()
967 sage: t = mu.parent().gen()
968 sage: mu + int(J.is_trivial())*(t-1)
971 The degree of an element is (by one definition) the degree
972 of its minimal polynomial::
974 sage: x = random_eja().random_element()
975 sage: x.degree() == x.minimal_polynomial().degree()
978 The minimal polynomial and the characteristic polynomial coincide
979 and are known (see Alizadeh, Example 11.11) for all elements of
980 the spin factor algebra that aren't scalar multiples of the
981 identity. We require the dimension of the algebra to be at least
982 two here so that said elements actually exist::
984 sage: d_max = JordanSpinEJA._max_random_instance_dimension()
985 sage: n = ZZ.random_element(2, max(2,d_max))
986 sage: J = JordanSpinEJA(n)
987 sage: y = J.random_element()
988 sage: while y == y.coefficient(0)*J.one():
989 ....: y = J.random_element()
990 sage: y0 = y.to_vector()[0]
991 sage: y_bar = y.to_vector()[1:]
992 sage: actual = y.minimal_polynomial()
993 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
994 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
995 sage: bool(actual == expected)
998 The minimal polynomial should always kill its element::
1000 sage: x = random_eja().random_element() # long time
1001 sage: p = x.minimal_polynomial() # long time
1002 sage: x.apply_univariate_polynomial(p) # long time
1005 The minimal polynomial is invariant under a change of basis,
1006 and in particular, a re-scaling of the basis::
1008 sage: d_max = RealSymmetricEJA._max_random_instance_dimension()
1009 sage: d = ZZ.random_element(1, d_max)
1010 sage: n = RealSymmetricEJA._max_random_instance_size(d)
1011 sage: J1 = RealSymmetricEJA(n)
1012 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
1013 sage: X = random_matrix(AA,n)
1014 sage: X = X*X.transpose()
1017 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
1022 # Pretty sure we know what the minimal polynomial of
1023 # the zero operator is going to be. This ensures
1024 # consistency of e.g. the polynomial variable returned
1025 # in the "normal" case without us having to think about it.
1026 return self
.operator().minimal_polynomial()
1028 # If we don't orthonormalize the subalgebra's basis, then the
1029 # first two monomials in the subalgebra will be self^0 and
1030 # self^1... assuming that self^1 is not a scalar multiple of
1031 # self^0 (the unit element). We special case these to avoid
1032 # having to solve a system to coerce self into the subalgebra.
1033 A
= self
.subalgebra_generated_by(orthonormalize
=False)
1035 if A
.dimension() == 1:
1036 # Does a solve to find the scalar multiple alpha such that
1037 # alpha*unit = self. We have to do this because the basis
1038 # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
1039 unit
= self
.parent().one()
1040 alpha
= self
.to_vector() / unit
.to_vector()
1041 return (unit
.operator()*alpha
).minimal_polynomial()
1043 # If the dimension of the subalgebra is >= 2, then we just
1044 # use the second basis element.
1045 return A
.monomial(1).operator().minimal_polynomial()
1049 def to_matrix(self
):
1051 Return an (often more natural) representation of this element as a
1054 Every finite-dimensional Euclidean Jordan Algebra is a direct
1055 sum of five simple algebras, four of which comprise Hermitian
1056 matrices. This method returns a "natural" matrix
1057 representation of this element as either a Hermitian matrix or
1062 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1064 ....: QuaternionHermitianEJA,
1065 ....: RealSymmetricEJA)
1069 sage: J = ComplexHermitianEJA(3)
1072 sage: J.one().to_matrix()
1083 sage: J = QuaternionHermitianEJA(2)
1086 sage: J.one().to_matrix()
1093 This also works in Cartesian product algebras::
1095 sage: J1 = HadamardEJA(1)
1096 sage: J2 = RealSymmetricEJA(2)
1097 sage: J = cartesian_product([J1,J2])
1098 sage: x = sum(J.gens())
1099 sage: x.to_matrix()[0]
1101 sage: x.to_matrix()[1]
1102 [ 1 0.7071067811865475?]
1103 [0.7071067811865475? 1]
1106 B
= self
.parent().matrix_basis()
1107 W
= self
.parent().matrix_space()
1109 # This is just a manual "from_vector()", but of course
1110 # matrix spaces aren't vector spaces in sage, so they
1111 # don't have a from_vector() method.
1112 return W
.linear_combination( zip(B
, self
.to_vector()) )
1118 The norm of this element with respect to :meth:`inner_product`.
1122 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1127 sage: J = HadamardEJA(2)
1128 sage: x = sum(J.gens())
1134 sage: J = JordanSpinEJA(4)
1135 sage: x = sum(J.gens())
1140 return self
.inner_product(self
).sqrt()
1145 Return the left-multiplication-by-this-element
1146 operator on the ambient algebra.
1150 sage: from mjo.eja.eja_algebra import random_eja
1154 sage: J = random_eja()
1155 sage: x,y = J.random_elements(2)
1156 sage: x.operator()(y) == x*y
1158 sage: y.operator()(x) == x*y
1163 left_mult_by_self
= lambda y
: self
*y
1164 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1165 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1168 def quadratic_representation(self
, other
=None):
1170 Return the quadratic representation of this element.
1174 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1179 The explicit form in the spin factor algebra is given by
1180 Alizadeh's Example 11.12::
1182 sage: x = JordanSpinEJA.random_instance().random_element()
1183 sage: x_vec = x.to_vector()
1184 sage: Q = matrix.identity(x.base_ring(), 0)
1185 sage: n = x_vec.degree()
1188 ....: x_bar = x_vec[1:]
1189 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1190 ....: B = 2*x0*x_bar.row()
1191 ....: C = 2*x0*x_bar.column()
1192 ....: D = matrix.identity(x.base_ring(), n-1)
1193 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1194 ....: D = D + 2*x_bar.tensor_product(x_bar)
1195 ....: Q = matrix.block(2,2,[A,B,C,D])
1196 sage: Q == x.quadratic_representation().matrix()
1199 Test all of the properties from Theorem 11.2 in Alizadeh::
1201 sage: J = random_eja()
1202 sage: x,y = J.random_elements(2)
1203 sage: Lx = x.operator()
1204 sage: Lxx = (x*x).operator()
1205 sage: Qx = x.quadratic_representation()
1206 sage: Qy = y.quadratic_representation()
1207 sage: Qxy = x.quadratic_representation(y)
1208 sage: Qex = J.one().quadratic_representation(x)
1209 sage: n = ZZ.random_element(10)
1210 sage: Qxn = (x^n).quadratic_representation()
1214 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1217 Property 2 (multiply on the right for :trac:`28272`):
1219 sage: alpha = J.base_ring().random_element()
1220 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1225 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1228 sage: not x.is_invertible() or (
1231 ....: x.inverse().quadratic_representation() )
1234 sage: Qxy(J.one()) == x*y
1239 sage: not x.is_invertible() or (
1240 ....: x.quadratic_representation(x.inverse())*Qx
1241 ....: == Qx*x.quadratic_representation(x.inverse()) )
1244 sage: not x.is_invertible() or (
1245 ....: x.quadratic_representation(x.inverse())*Qx
1247 ....: 2*Lx*Qex - Qx )
1250 sage: 2*Lx*Qex - Qx == Lxx
1255 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1265 sage: not x.is_invertible() or (
1266 ....: Qx*x.inverse().operator() == Lx )
1271 sage: not x.operator_commutes_with(y) or (
1272 ....: Qx(y)^n == Qxn(y^n) )
1278 elif not other
in self
.parent():
1279 raise TypeError("'other' must live in the same algebra")
1282 M
= other
.operator()
1283 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1287 def spectral_decomposition(self
):
1289 Return the unique spectral decomposition of this element.
1293 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1294 element's left-multiplication-by operator to the subalgebra it
1295 generates. We then compute the spectral decomposition of that
1296 operator, and the spectral projectors we get back must be the
1297 left-multiplication-by operators for the idempotents we
1298 seek. Thus applying them to the identity element gives us those
1301 Since the eigenvalues are required to be distinct, we take
1302 the spectral decomposition of the zero element to be zero
1303 times the identity element of the algebra (which is idempotent,
1308 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1312 The spectral decomposition of the identity is ``1`` times itself,
1313 and the spectral decomposition of zero is ``0`` times the identity::
1315 sage: J = RealSymmetricEJA(3)
1318 sage: J.one().spectral_decomposition()
1320 sage: J.zero().spectral_decomposition()
1325 sage: J = RealSymmetricEJA(4)
1326 sage: x = sum(J.gens())
1327 sage: sd = x.spectral_decomposition()
1332 sage: c0.inner_product(c1) == 0
1334 sage: c0.is_idempotent()
1336 sage: c1.is_idempotent()
1338 sage: c0 + c1 == J.one()
1340 sage: l0*c0 + l1*c1 == x
1343 The spectral decomposition should work in subalgebras, too::
1345 sage: J = RealSymmetricEJA(4)
1346 sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
1347 sage: A = 2*b5 - 2*b8
1348 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1349 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1350 sage: (f0, f1, f2) = J1.gens()
1351 sage: f0.spectral_decomposition()
1352 [(0, 1.000000000000000?*c2), (1, 1.000000000000000?*c0)]
1355 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1357 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1358 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1361 def subalgebra_generated_by(self
, **kwargs
):
1363 Return the associative subalgebra of the parent EJA generated
1366 Since our parent algebra is unital, we want "subalgebra" to mean
1367 "unital subalgebra" as well; thus the subalgebra that an element
1368 generates will itself be a Euclidean Jordan algebra after
1369 restricting the algebra operations appropriately. This is the
1370 subalgebra that Faraut and Korányi work with in section II.2, for
1375 sage: from mjo.eja.eja_algebra import (random_eja,
1377 ....: RealSymmetricEJA)
1381 We can create subalgebras of Cartesian product EJAs that are not
1382 themselves Cartesian product EJAs (they're just "regular" EJAs)::
1384 sage: J1 = HadamardEJA(3)
1385 sage: J2 = RealSymmetricEJA(2)
1386 sage: J = cartesian_product([J1,J2])
1387 sage: J.one().subalgebra_generated_by()
1388 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
1392 This subalgebra, being composed of only powers, is associative::
1394 sage: x0 = random_eja().random_element()
1395 sage: A = x0.subalgebra_generated_by(orthonormalize=False)
1396 sage: x,y,z = A.random_elements(3)
1397 sage: (x*y)*z == x*(y*z)
1400 Squaring in the subalgebra should work the same as in
1403 sage: x = random_eja().random_element()
1404 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1405 sage: A(x^2) == A(x)*A(x)
1408 By definition, the subalgebra generated by the zero element is
1409 the one-dimensional algebra generated by the identity
1410 element... unless the original algebra was trivial, in which
1411 case the subalgebra is trivial too::
1413 sage: A = random_eja().zero().subalgebra_generated_by()
1414 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1418 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1419 A
= self
.parent().subalgebra(powers
,
1424 A
.one
.set_cache(A(self
.parent().one()))
1428 def subalgebra_idempotent(self
):
1430 Find an idempotent in the associative subalgebra I generate
1431 using Proposition 2.3.5 in Baes.
1435 sage: from mjo.eja.eja_algebra import random_eja
1439 Ensure that we can find an idempotent in a non-trivial algebra
1440 where there are non-nilpotent elements, or that we get the dumb
1441 solution in the trivial algebra::
1443 sage: J = random_eja(field=QQ, orthonormalize=False)
1444 sage: x = J.random_element()
1445 sage: while x.is_nilpotent() and not J.is_trivial():
1446 ....: x = J.random_element()
1447 sage: c = x.subalgebra_idempotent()
1452 if self
.parent().is_trivial():
1455 if self
.is_nilpotent():
1456 raise ValueError("this only works with non-nilpotent elements!")
1458 J
= self
.subalgebra_generated_by()
1461 # The image of the matrix of left-u^m-multiplication
1462 # will be minimal for some natural number s...
1464 minimal_dim
= J
.dimension()
1465 for i
in range(1, minimal_dim
):
1466 this_dim
= (u
**i
).operator().matrix().image().dimension()
1467 if this_dim
< minimal_dim
:
1468 minimal_dim
= this_dim
1471 # Now minimal_matrix should correspond to the smallest
1472 # non-zero subspace in Baes's (or really, Koecher's)
1475 # However, we need to restrict the matrix to work on the
1476 # subspace... or do we? Can't we just solve, knowing that
1477 # A(c) = u^(s+1) should have a solution in the big space,
1480 # Beware, solve_right() means that we're using COLUMN vectors.
1481 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1483 A
= u_next
.operator().matrix()
1484 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1486 # Now c is the idempotent we want, but it still lives in the subalgebra.
1487 return c
.superalgebra_element()
1492 Return my trace, the sum of my eigenvalues.
1494 In a trivial algebra, however you want to look at it, the trace is
1495 an empty sum for which we declare the result to be zero.
1499 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1506 sage: J = TrivialEJA()
1507 sage: J.zero().trace()
1511 sage: J = JordanSpinEJA(3)
1512 sage: x = sum(J.gens())
1518 sage: J = HadamardEJA(5)
1519 sage: J.one().trace()
1524 The trace of an element is a real number::
1526 sage: J = random_eja()
1527 sage: J.random_element().trace() in RLF
1530 The trace is linear::
1532 sage: J = random_eja()
1533 sage: x,y = J.random_elements(2)
1534 sage: alpha = J.base_ring().random_element()
1535 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1538 The trace of a square is nonnegative::
1540 sage: x = random_eja().random_element()
1541 sage: (x*x).trace() >= 0
1549 # Special case for the trivial algebra where
1550 # the trace is an empty sum.
1551 return P
.base_ring().zero()
1553 p
= P
._charpoly
_coefficients
()[r
-1]
1554 # The _charpoly_coeff function already adds the factor of
1555 # -1 to ensure that _charpoly_coeff(r-1) is really what
1556 # appears in front of t^{r-1} in the charpoly. However,
1557 # we want the negative of THAT for the trace.
1558 return -p(*self
.to_vector())
1560 def operator_inner_product(self
, other
):
1562 Return the operator inner product of myself and ``other``.
1564 The "operator inner product," whose name is not standard, is
1565 defined be the usual linear-algebraic trace of the
1566 ``(x*y).operator()``.
1568 Proposition III.1.5 in Faraut and Korányi shows that on any
1569 Euclidean Jordan algebra, this is another associative inner
1570 product under which the cone of squares is symmetric.
1572 This *probably* works even if the basis hasn't been
1573 orthonormalized because the eigenvalues of the corresponding
1574 matrix don't change when the basis does (they're preserved by
1575 any similarity transformation).
1579 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1580 ....: RealSymmetricEJA,
1581 ....: ComplexHermitianEJA,
1586 Proposition III.4.2 of Faraut and Korányi shows that on a
1587 simple algebra of rank `r` and dimension `n`, this inner
1588 product is `n/r` times the canonical
1589 :meth:`trace_inner_product`::
1591 sage: J = JordanSpinEJA(4, field=QQ)
1592 sage: x,y = J.random_elements(2)
1593 sage: n = J.dimension()
1595 sage: actual = x.operator_inner_product(y)
1596 sage: expected = (n/r)*x.trace_inner_product(y)
1597 sage: actual == expected
1602 sage: J = RealSymmetricEJA(3)
1603 sage: x,y = J.random_elements(2)
1604 sage: n = J.dimension()
1606 sage: actual = x.operator_inner_product(y)
1607 sage: expected = (n/r)*x.trace_inner_product(y)
1608 sage: actual == expected
1613 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
1614 sage: x,y = J.random_elements(2)
1615 sage: n = J.dimension()
1617 sage: actual = x.operator_inner_product(y)
1618 sage: expected = (n/r)*x.trace_inner_product(y)
1619 sage: actual == expected
1624 The operator inner product is commutative, bilinear, and
1627 sage: J = random_eja()
1628 sage: x,y,z = J.random_elements(3)
1630 sage: x.operator_inner_product(y) == y.operator_inner_product(x)
1633 sage: a = J.base_ring().random_element()
1634 sage: actual = (a*(x+z)).operator_inner_product(y)
1635 sage: expected = ( a*x.operator_inner_product(y) +
1636 ....: a*z.operator_inner_product(y) )
1637 sage: actual == expected
1639 sage: actual = x.operator_inner_product(a*(y+z))
1640 sage: expected = ( a*x.operator_inner_product(y) +
1641 ....: a*x.operator_inner_product(z) )
1642 sage: actual == expected
1645 sage: actual = (x*y).operator_inner_product(z)
1646 sage: expected = y.operator_inner_product(x*z)
1647 sage: actual == expected
1651 if not other
in self
.parent():
1652 raise TypeError("'other' must live in the same algebra")
1654 return (self
*other
).operator().matrix().trace()
1657 def trace_inner_product(self
, other
):
1659 Return the trace inner product of myself and ``other``.
1663 sage: from mjo.eja.eja_algebra import random_eja
1667 The trace inner product is commutative, bilinear, and associative::
1669 sage: J = random_eja()
1670 sage: x,y,z = J.random_elements(3)
1672 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1675 sage: a = J.base_ring().random_element()
1676 sage: actual = (a*(x+z)).trace_inner_product(y)
1677 sage: expected = ( a*x.trace_inner_product(y) +
1678 ....: a*z.trace_inner_product(y) )
1679 sage: actual == expected
1681 sage: actual = x.trace_inner_product(a*(y+z))
1682 sage: expected = ( a*x.trace_inner_product(y) +
1683 ....: a*x.trace_inner_product(z) )
1684 sage: actual == expected
1687 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1691 if not other
in self
.parent():
1692 raise TypeError("'other' must live in the same algebra")
1694 return (self
*other
).trace()
1697 def trace_norm(self
):
1699 The norm of this element with respect to :meth:`trace_inner_product`.
1703 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1708 sage: J = HadamardEJA(2)
1709 sage: x = sum(J.gens())
1710 sage: x.trace_norm()
1715 sage: J = JordanSpinEJA(4)
1716 sage: x = sum(J.gens())
1717 sage: x.trace_norm()
1721 return self
.trace_inner_product(self
).sqrt()
1724 class CartesianProductEJAElement(FiniteDimensionalEJAElement
):
1727 Compute the determinant of this product-element using the
1728 determianants of its factors.
1730 This result Follows from the spectral decomposition of (say)
1731 the pair `(x,y)` in terms of the Jordan frame `\left\{ (c_1,
1732 0),(c_2, 0),...,(0,d_1),(0,d_2),... \right\}.
1734 from sage
.misc
.misc_c
import prod
1735 return prod( f
.det() for f
in self
.cartesian_factors() )
1737 def to_matrix(self
):
1738 # An override is necessary to call our custom _scale().
1739 B
= self
.parent().matrix_basis()
1740 W
= self
.parent().matrix_space()
1742 # Aaaaand linear combinations don't work in Cartesian
1743 # product spaces, even though they provide a method with
1744 # that name. This is hidden behind an "if" because the
1745 # _scale() function is slow.
1746 pairs
= zip(B
, self
.to_vector())
1747 return W
.sum( _scale(b
, alpha
) for (b
,alpha
) in pairs
)