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)
293 sage: Lx = x.operator()
294 sage: Ly = y.operator()
295 sage: Lz = z.operator()
296 sage: Lzy = (z*y).operator()
297 sage: Lxy = (x*y).operator()
298 sage: Lxz = (x*z).operator()
299 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
302 Test the third polarization identity from my notes or from
305 sage: u,y,z = random_eja().random_elements(3)
306 sage: Lu = u.operator()
307 sage: Ly = y.operator()
308 sage: Lz = z.operator()
309 sage: Lzy = (z*y).operator()
310 sage: Luy = (u*y).operator()
311 sage: Luz = (u*z).operator()
312 sage: Luyz = (u*(y*z)).operator()
313 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
314 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
315 sage: bool(lhs == rhs)
319 if not other
in self
.parent():
320 raise TypeError("'other' must live in the same algebra")
329 Return my determinant, the product of my eigenvalues.
333 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
335 ....: RealSymmetricEJA,
336 ....: ComplexHermitianEJA,
341 sage: J = JordanSpinEJA(2)
342 sage: x = sum( J.gens() )
348 sage: J = JordanSpinEJA(3)
349 sage: x = sum( J.gens() )
353 The determinant of the sole element in the rank-zero trivial
354 algebra is ``1``, by three paths of reasoning. First, its
355 characteristic polynomial is a constant ``1``, so the constant
356 term in that polynomial is ``1``. Second, the characteristic
357 polynomial evaluated at zero is again ``1``. And finally, the
358 (empty) product of its eigenvalues is likewise just unity::
360 sage: J = TrivialEJA()
366 An element is invertible if and only if its determinant is
369 sage: x = random_eja().random_element()
370 sage: x.is_invertible() == (x.det() != 0)
373 Ensure that the determinant is multiplicative on an associative
374 subalgebra as in Faraut and Korányi's Proposition II.2.2::
376 sage: J = random_eja().random_element().subalgebra_generated_by()
377 sage: x,y = J.random_elements(2)
378 sage: (x*y).det() == x.det()*y.det()
381 The determinant in real matrix algebras is the usual determinant::
383 sage: X = matrix.random(QQ,3)
385 sage: J1 = RealSymmetricEJA(3)
386 sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
387 sage: expected = X.det()
388 sage: actual1 = J1(X).det()
389 sage: actual2 = J2(X).det()
390 sage: actual1 == expected
392 sage: actual2 == expected
400 # Special case, since we don't get the a0=1
401 # coefficient when the rank of the algebra
403 return P
.base_ring().one()
405 p
= P
._charpoly
_coefficients
()[0]
406 # The _charpoly_coeff function already adds the factor of -1
407 # to ensure that _charpoly_coefficients()[0] is really what
408 # appears in front of t^{0} in the charpoly. However, we want
409 # (-1)^r times THAT for the determinant.
410 return ((-1)**r
)*p(*self
.to_vector())
416 Return the Jordan-multiplicative inverse of this element.
420 In general we appeal to the quadratic representation as in
421 Koecher's Theorem 12 in Chapter III, Section 5. But if the
422 parent algebra's "characteristic polynomial of" coefficients
423 happen to be cached, then we use Proposition II.2.4 in Faraut
424 and Korányi which gives a formula for the inverse based on the
425 characteristic polynomial and the Cayley-Hamilton theorem for
426 Euclidean Jordan algebras::
430 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
436 The inverse in the spin factor algebra is given in Alizadeh's
439 sage: J = JordanSpinEJA.random_instance()
440 sage: x = J.random_element()
441 sage: while not x.is_invertible():
442 ....: x = J.random_element()
443 sage: x_vec = x.to_vector()
445 sage: x_bar = x_vec[1:]
446 sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
447 sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
448 sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
449 sage: x.inverse() == J.from_vector(x_inverse)
452 Trying to invert a non-invertible element throws an error:
454 sage: JordanSpinEJA(3).zero().inverse()
455 Traceback (most recent call last):
457 ZeroDivisionError: element is not invertible
461 The identity element is its own inverse::
463 sage: J = random_eja()
464 sage: J.one().inverse() == J.one()
467 If an element has an inverse, it acts like one::
469 sage: J = random_eja()
470 sage: x = J.random_element()
471 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
474 The inverse of the inverse is what we started with::
476 sage: J = random_eja()
477 sage: x = J.random_element()
478 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
481 Proposition II.2.3 in Faraut and Korányi says that the inverse
482 of an element is the inverse of its left-multiplication operator
483 applied to the algebra's identity, when that inverse exists::
485 sage: J = random_eja()
486 sage: x = J.random_element()
487 sage: (not x.operator().is_invertible()) or (
488 ....: x.operator().inverse()(J.one()) == x.inverse() )
491 Check that the fast (cached) and slow algorithms give the same
494 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
495 sage: x = J.random_element() # long time
496 sage: while not x.is_invertible(): # long time
497 ....: x = J.random_element() # long time
498 sage: slow = x.inverse() # long time
499 sage: _ = J._charpoly_coefficients() # long time
500 sage: fast = x.inverse() # long time
501 sage: slow == fast # long time
504 not_invertible_msg
= "element is not invertible"
505 if self
.parent()._charpoly
_coefficients
.is_in_cache():
506 # We can invert using our charpoly if it will be fast to
507 # compute. If the coefficients are cached, our rank had
509 if self
.det().is_zero():
510 raise ZeroDivisionError(not_invertible_msg
)
511 r
= self
.parent().rank()
512 a
= self
.characteristic_polynomial().coefficients(sparse
=False)
513 return (-1)**(r
+1)*sum(a
[i
+1]*self
**i
for i
in range(r
))/self
.det()
516 inv
= (~self
.quadratic_representation())(self
)
517 self
.is_invertible
.set_cache(True)
519 except ZeroDivisionError:
520 self
.is_invertible
.set_cache(False)
521 raise ZeroDivisionError(not_invertible_msg
)
525 def is_invertible(self
):
527 Return whether or not this element is invertible.
531 If computing my determinant will be fast, we do so and compare
532 with zero (Proposition II.2.4 in Faraut and
533 Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
534 reduces the problem to the invertibility of my quadratic
539 sage: from mjo.eja.eja_algebra import random_eja
543 The identity element is always invertible::
545 sage: J = random_eja()
546 sage: J.one().is_invertible()
549 The zero element is never invertible in a non-trivial algebra::
551 sage: J = random_eja()
552 sage: (not J.is_trivial()) and J.zero().is_invertible()
555 Test that the fast (cached) and slow algorithms give the same
558 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
559 sage: x = J.random_element() # long time
560 sage: slow = x.is_invertible() # long time
561 sage: _ = J._charpoly_coefficients() # long time
562 sage: fast = x.is_invertible() # long time
563 sage: slow == fast # long time
567 if self
.parent().is_trivial():
572 if self
.parent()._charpoly
_coefficients
.is_in_cache():
573 # The determinant will be quicker than inverting the
574 # quadratic representation, most likely.
575 return (not self
.det().is_zero())
577 # The easiest way to determine if I'm invertible is to try.
579 inv
= (~self
.quadratic_representation())(self
)
580 self
.inverse
.set_cache(inv
)
582 except ZeroDivisionError:
586 def is_primitive_idempotent(self
):
588 Return whether or not this element is a primitive (or minimal)
591 A primitive idempotent is a non-zero idempotent that is not
592 the sum of two other non-zero idempotents. Remark 2.7.15 in
593 Baes shows that this is what he refers to as a "minimal
596 An element of a Euclidean Jordan algebra is a minimal idempotent
597 if it :meth:`is_idempotent` and if its Peirce subalgebra
598 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
603 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
604 ....: RealSymmetricEJA,
610 This method is sloooooow.
614 The spectral decomposition of a non-regular element should always
615 contain at least one non-minimal idempotent::
617 sage: J = RealSymmetricEJA(3)
618 sage: x = sum(J.gens())
621 sage: [ c.is_primitive_idempotent()
622 ....: for (l,c) in x.spectral_decomposition() ]
625 On the other hand, the spectral decomposition of a regular
626 element should always be in terms of minimal idempotents::
628 sage: J = JordanSpinEJA(4)
629 sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
632 sage: [ c.is_primitive_idempotent()
633 ....: for (l,c) in x.spectral_decomposition() ]
638 The identity element is minimal only in an EJA of rank one::
640 sage: J = random_eja()
641 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
644 A non-idempotent cannot be a minimal idempotent::
646 sage: J = JordanSpinEJA(4)
647 sage: x = J.random_element()
648 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
651 Proposition 2.7.19 in Baes says that an element is a minimal
652 idempotent if and only if it's idempotent with trace equal to
655 sage: J = JordanSpinEJA(4)
656 sage: x = J.random_element()
657 sage: expected = (x.is_idempotent() and x.trace() == 1)
658 sage: actual = x.is_primitive_idempotent()
659 sage: actual == expected
662 Primitive idempotents must be non-zero::
664 sage: J = random_eja()
665 sage: J.zero().is_idempotent()
667 sage: J.zero().is_primitive_idempotent()
670 As a consequence of the fact that primitive idempotents must
671 be non-zero, there are no primitive idempotents in a trivial
672 Euclidean Jordan algebra::
674 sage: J = TrivialEJA()
675 sage: J.one().is_idempotent()
677 sage: J.one().is_primitive_idempotent()
681 if not self
.is_idempotent():
687 (_
,_
,J1
) = self
.parent().peirce_decomposition(self
)
688 return (J1
.dimension() == 1)
691 def is_nilpotent(self
):
693 Return whether or not some power of this element is zero.
697 We use Theorem 5 in Chapter III of Koecher, which says that
698 an element ``x`` is nilpotent if and only if ``x.operator()``
699 is nilpotent. And it is a basic fact of linear algebra that
700 an operator on an `n`-dimensional space is nilpotent if and
701 only if, when raised to the `n`th power, it equals the zero
702 operator (for example, see Axler Corollary 8.8).
706 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
711 sage: J = JordanSpinEJA(3)
712 sage: x = sum(J.gens())
713 sage: x.is_nilpotent()
718 The identity element is never nilpotent, except in a trivial EJA::
720 sage: J = random_eja()
721 sage: J.one().is_nilpotent() and not J.is_trivial()
724 The additive identity is always nilpotent::
726 sage: random_eja().zero().is_nilpotent()
731 zero_operator
= P
.zero().operator()
732 return self
.operator()**P
.dimension() == zero_operator
735 def is_regular(self
):
737 Return whether or not this is a regular element.
741 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
746 The identity element always has degree one, but any element
747 linearly-independent from it is regular::
749 sage: J = JordanSpinEJA(5)
750 sage: J.one().is_regular()
752 sage: b0, b1, b2, b3, b4 = J.gens()
755 sage: for x in J.gens():
756 ....: (J.one() + x).is_regular()
765 The zero element should never be regular, unless the parent
766 algebra has dimension less than or equal to one::
768 sage: J = random_eja()
769 sage: J.dimension() <= 1 or not J.zero().is_regular()
772 The unit element isn't regular unless the algebra happens to
773 consist of only its scalar multiples::
775 sage: J = random_eja()
776 sage: J.dimension() <= 1 or not J.one().is_regular()
780 return self
.degree() == self
.parent().rank()
785 Return the degree of this element, which is defined to be
786 the degree of its minimal polynomial.
794 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
799 sage: J = JordanSpinEJA(4)
800 sage: J.one().degree()
802 sage: b0,b1,b2,b3 = J.gens()
803 sage: (b0 - b1).degree()
806 In the spin factor algebra (of rank two), all elements that
807 aren't multiples of the identity are regular::
809 sage: J = JordanSpinEJA.random_instance()
810 sage: n = J.dimension()
811 sage: x = J.random_element()
812 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
817 The zero and unit elements are both of degree one in nontrivial
820 sage: J = random_eja()
821 sage: d = J.zero().degree()
822 sage: (J.is_trivial() and d == 0) or d == 1
824 sage: d = J.one().degree()
825 sage: (J.is_trivial() and d == 0) or d == 1
828 Our implementation agrees with the definition::
830 sage: x = random_eja().random_element()
831 sage: x.degree() == x.minimal_polynomial().degree()
835 n
= self
.parent().dimension()
838 # The minimal polynomial is an empty product, i.e. the
839 # constant polynomial "1" having degree zero.
842 # The minimal polynomial of zero in a nontrivial algebra
843 # is "t", and is of degree one.
846 # If this is a nonzero element of a nontrivial algebra, it
847 # has degree at least one. It follows that, in an algebra
848 # of dimension one, the degree must be actually one.
851 # BEWARE: The subalgebra_generated_by() method uses the result
852 # of this method to construct a basis for the subalgebra. That
853 # means, in particular, that we cannot implement this method
854 # as ``self.subalgebra_generated_by().dimension()``.
856 # Algorithm: keep appending (vector representations of) powers
857 # self as rows to a matrix and echelonizing it. When its rank
858 # stops increasing, we've reached a redundancy.
860 # Given the special cases above, we can assume that "self" is
861 # nonzero, the algebra is nontrivial, and that its dimension
863 M
= matrix([(self
.parent().one()).to_vector()])
866 # Specifying the row-reduction algorithm can e.g. help over
867 # AA because it avoids the RecursionError that gets thrown
868 # when we have to look too hard for a root.
870 # Beware: QQ supports an entirely different set of "algorithm"
871 # keywords than do AA and RR.
873 from sage
.rings
.all
import QQ
874 if self
.parent().base_ring() is not QQ
:
875 algo
= "scaled_partial_pivoting"
878 M
= matrix(M
.rows() + [(self
**d
).to_vector()])
881 if new_rank
== old_rank
:
890 def left_matrix(self
):
892 Our parent class defines ``left_matrix`` and ``matrix``
893 methods whose names are misleading. We don't want them.
895 raise NotImplementedError("use operator().matrix() instead")
900 def minimal_polynomial(self
):
902 Return the minimal polynomial of this element,
903 as a function of the variable `t`.
907 We restrict ourselves to the associative subalgebra
908 generated by this element, and then return the minimal
909 polynomial of this element's operator matrix (in that
910 subalgebra). This works by Baes Proposition 2.3.16.
914 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
915 ....: RealSymmetricEJA,
921 Keeping in mind that the polynomial ``1`` evaluates the identity
922 element (also the zero element) of the trivial algebra, it is clear
923 that the polynomial ``1`` is the minimal polynomial of the only
924 element in a trivial algebra::
926 sage: J = TrivialEJA()
927 sage: J.one().minimal_polynomial()
929 sage: J.zero().minimal_polynomial()
934 The minimal polynomial of the identity and zero elements are
935 always the same, except in trivial algebras where the minimal
936 polynomial of the unit/zero element is ``1``::
938 sage: J = random_eja()
939 sage: mu = J.one().minimal_polynomial()
940 sage: t = mu.parent().gen()
941 sage: mu + int(J.is_trivial())*(t-2)
943 sage: mu = J.zero().minimal_polynomial()
944 sage: t = mu.parent().gen()
945 sage: mu + int(J.is_trivial())*(t-1)
948 The degree of an element is (by one definition) the degree
949 of its minimal polynomial::
951 sage: x = random_eja().random_element()
952 sage: x.degree() == x.minimal_polynomial().degree()
955 The minimal polynomial and the characteristic polynomial coincide
956 and are known (see Alizadeh, Example 11.11) for all elements of
957 the spin factor algebra that aren't scalar multiples of the
958 identity. We require the dimension of the algebra to be at least
959 two here so that said elements actually exist::
961 sage: d_max = JordanSpinEJA._max_random_instance_dimension()
962 sage: n = ZZ.random_element(2, max(2,d_max))
963 sage: J = JordanSpinEJA(n)
964 sage: y = J.random_element()
965 sage: while y == y.coefficient(0)*J.one():
966 ....: y = J.random_element()
967 sage: y0 = y.to_vector()[0]
968 sage: y_bar = y.to_vector()[1:]
969 sage: actual = y.minimal_polynomial()
970 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
971 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
972 sage: bool(actual == expected)
975 The minimal polynomial should always kill its element::
977 sage: x = random_eja().random_element()
978 sage: p = x.minimal_polynomial()
979 sage: x.apply_univariate_polynomial(p)
982 The minimal polynomial is invariant under a change of basis,
983 and in particular, a re-scaling of the basis::
985 sage: d_max = RealSymmetricEJA._max_random_instance_dimension()
986 sage: d = ZZ.random_element(1, d_max)
987 sage: n = RealSymmetricEJA._max_random_instance_size(d)
988 sage: J1 = RealSymmetricEJA(n)
989 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
990 sage: X = random_matrix(AA,n)
991 sage: X = X*X.transpose()
994 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
999 # Pretty sure we know what the minimal polynomial of
1000 # the zero operator is going to be. This ensures
1001 # consistency of e.g. the polynomial variable returned
1002 # in the "normal" case without us having to think about it.
1003 return self
.operator().minimal_polynomial()
1005 # If we don't orthonormalize the subalgebra's basis, then the
1006 # first two monomials in the subalgebra will be self^0 and
1007 # self^1... assuming that self^1 is not a scalar multiple of
1008 # self^0 (the unit element). We special case these to avoid
1009 # having to solve a system to coerce self into the subalgebra.
1010 A
= self
.subalgebra_generated_by(orthonormalize
=False)
1012 if A
.dimension() == 1:
1013 # Does a solve to find the scalar multiple alpha such that
1014 # alpha*unit = self. We have to do this because the basis
1015 # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
1016 unit
= self
.parent().one()
1017 alpha
= self
.to_vector() / unit
.to_vector()
1018 return (unit
.operator()*alpha
).minimal_polynomial()
1020 # If the dimension of the subalgebra is >= 2, then we just
1021 # use the second basis element.
1022 return A
.monomial(1).operator().minimal_polynomial()
1026 def to_matrix(self
):
1028 Return an (often more natural) representation of this element as a
1031 Every finite-dimensional Euclidean Jordan Algebra is a direct
1032 sum of five simple algebras, four of which comprise Hermitian
1033 matrices. This method returns a "natural" matrix
1034 representation of this element as either a Hermitian matrix or
1039 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1041 ....: QuaternionHermitianEJA,
1042 ....: RealSymmetricEJA)
1046 sage: J = ComplexHermitianEJA(3)
1049 sage: J.one().to_matrix()
1060 sage: J = QuaternionHermitianEJA(2)
1063 sage: J.one().to_matrix()
1070 This also works in Cartesian product algebras::
1072 sage: J1 = HadamardEJA(1)
1073 sage: J2 = RealSymmetricEJA(2)
1074 sage: J = cartesian_product([J1,J2])
1075 sage: x = sum(J.gens())
1076 sage: x.to_matrix()[0]
1078 sage: x.to_matrix()[1]
1079 [ 1 0.7071067811865475?]
1080 [0.7071067811865475? 1]
1083 B
= self
.parent().matrix_basis()
1084 W
= self
.parent().matrix_space()
1086 # This is just a manual "from_vector()", but of course
1087 # matrix spaces aren't vector spaces in sage, so they
1088 # don't have a from_vector() method.
1089 return W
.linear_combination( zip(B
, self
.to_vector()) )
1095 The norm of this element with respect to :meth:`inner_product`.
1099 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1104 sage: J = HadamardEJA(2)
1105 sage: x = sum(J.gens())
1111 sage: J = JordanSpinEJA(4)
1112 sage: x = sum(J.gens())
1117 return self
.inner_product(self
).sqrt()
1122 Return the left-multiplication-by-this-element
1123 operator on the ambient algebra.
1127 sage: from mjo.eja.eja_algebra import random_eja
1131 sage: J = random_eja()
1132 sage: x,y = J.random_elements(2)
1133 sage: x.operator()(y) == x*y
1135 sage: y.operator()(x) == x*y
1140 left_mult_by_self
= lambda y
: self
*y
1141 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1142 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1145 def quadratic_representation(self
, other
=None):
1147 Return the quadratic representation of this element.
1151 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1156 The explicit form in the spin factor algebra is given by
1157 Alizadeh's Example 11.12::
1159 sage: x = JordanSpinEJA.random_instance().random_element()
1160 sage: x_vec = x.to_vector()
1161 sage: Q = matrix.identity(x.base_ring(), 0)
1162 sage: n = x_vec.degree()
1165 ....: x_bar = x_vec[1:]
1166 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1167 ....: B = 2*x0*x_bar.row()
1168 ....: C = 2*x0*x_bar.column()
1169 ....: D = matrix.identity(x.base_ring(), n-1)
1170 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1171 ....: D = D + 2*x_bar.tensor_product(x_bar)
1172 ....: Q = matrix.block(2,2,[A,B,C,D])
1173 sage: Q == x.quadratic_representation().matrix()
1176 Test all of the properties from Theorem 11.2 in Alizadeh::
1178 sage: J = random_eja()
1179 sage: x,y = J.random_elements(2)
1180 sage: Lx = x.operator()
1181 sage: Lxx = (x*x).operator()
1182 sage: Qx = x.quadratic_representation()
1183 sage: Qy = y.quadratic_representation()
1184 sage: Qxy = x.quadratic_representation(y)
1185 sage: Qex = J.one().quadratic_representation(x)
1186 sage: n = ZZ.random_element(10)
1187 sage: Qxn = (x^n).quadratic_representation()
1191 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1194 Property 2 (multiply on the right for :trac:`28272`):
1196 sage: alpha = J.base_ring().random_element()
1197 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1202 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1205 sage: not x.is_invertible() or (
1208 ....: x.inverse().quadratic_representation() )
1211 sage: Qxy(J.one()) == x*y
1216 sage: not x.is_invertible() or (
1217 ....: x.quadratic_representation(x.inverse())*Qx
1218 ....: == Qx*x.quadratic_representation(x.inverse()) )
1221 sage: not x.is_invertible() or (
1222 ....: x.quadratic_representation(x.inverse())*Qx
1224 ....: 2*Lx*Qex - Qx )
1227 sage: 2*Lx*Qex - Qx == Lxx
1232 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1242 sage: not x.is_invertible() or (
1243 ....: Qx*x.inverse().operator() == Lx )
1248 sage: not x.operator_commutes_with(y) or (
1249 ....: Qx(y)^n == Qxn(y^n) )
1255 elif not other
in self
.parent():
1256 raise TypeError("'other' must live in the same algebra")
1259 M
= other
.operator()
1260 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1264 def spectral_decomposition(self
):
1266 Return the unique spectral decomposition of this element.
1270 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1271 element's left-multiplication-by operator to the subalgebra it
1272 generates. We then compute the spectral decomposition of that
1273 operator, and the spectral projectors we get back must be the
1274 left-multiplication-by operators for the idempotents we
1275 seek. Thus applying them to the identity element gives us those
1278 Since the eigenvalues are required to be distinct, we take
1279 the spectral decomposition of the zero element to be zero
1280 times the identity element of the algebra (which is idempotent,
1285 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1289 The spectral decomposition of the identity is ``1`` times itself,
1290 and the spectral decomposition of zero is ``0`` times the identity::
1292 sage: J = RealSymmetricEJA(3)
1295 sage: J.one().spectral_decomposition()
1297 sage: J.zero().spectral_decomposition()
1302 sage: J = RealSymmetricEJA(4)
1303 sage: x = sum(J.gens())
1304 sage: sd = x.spectral_decomposition()
1309 sage: c0.inner_product(c1) == 0
1311 sage: c0.is_idempotent()
1313 sage: c1.is_idempotent()
1315 sage: c0 + c1 == J.one()
1317 sage: l0*c0 + l1*c1 == x
1320 The spectral decomposition should work in subalgebras, too::
1322 sage: J = RealSymmetricEJA(4)
1323 sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
1324 sage: A = 2*b5 - 2*b8
1325 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1326 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1327 sage: (f0, f1, f2) = J1.gens()
1328 sage: f0.spectral_decomposition()
1329 [(0, 1.000000000000000?*c2), (1, 1.000000000000000?*c0)]
1332 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1334 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1335 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1338 def subalgebra_generated_by(self
, **kwargs
):
1340 Return the associative subalgebra of the parent EJA generated
1343 Since our parent algebra is unital, we want "subalgebra" to mean
1344 "unital subalgebra" as well; thus the subalgebra that an element
1345 generates will itself be a Euclidean Jordan algebra after
1346 restricting the algebra operations appropriately. This is the
1347 subalgebra that Faraut and Korányi work with in section II.2, for
1352 sage: from mjo.eja.eja_algebra import (random_eja,
1354 ....: RealSymmetricEJA)
1358 We can create subalgebras of Cartesian product EJAs that are not
1359 themselves Cartesian product EJAs (they're just "regular" EJAs)::
1361 sage: J1 = HadamardEJA(3)
1362 sage: J2 = RealSymmetricEJA(2)
1363 sage: J = cartesian_product([J1,J2])
1364 sage: J.one().subalgebra_generated_by()
1365 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
1369 This subalgebra, being composed of only powers, is associative::
1371 sage: x0 = random_eja().random_element()
1372 sage: A = x0.subalgebra_generated_by()
1373 sage: x,y,z = A.random_elements(3)
1374 sage: (x*y)*z == x*(y*z)
1377 Squaring in the subalgebra should work the same as in
1380 sage: x = random_eja().random_element()
1381 sage: A = x.subalgebra_generated_by()
1382 sage: A(x^2) == A(x)*A(x)
1385 By definition, the subalgebra generated by the zero element is
1386 the one-dimensional algebra generated by the identity
1387 element... unless the original algebra was trivial, in which
1388 case the subalgebra is trivial too::
1390 sage: A = random_eja().zero().subalgebra_generated_by()
1391 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1395 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1396 A
= self
.parent().subalgebra(powers
,
1401 A
.one
.set_cache(A(self
.parent().one()))
1405 def subalgebra_idempotent(self
):
1407 Find an idempotent in the associative subalgebra I generate
1408 using Proposition 2.3.5 in Baes.
1412 sage: from mjo.eja.eja_algebra import random_eja
1416 Ensure that we can find an idempotent in a non-trivial algebra
1417 where there are non-nilpotent elements, or that we get the dumb
1418 solution in the trivial algebra::
1420 sage: J = random_eja()
1421 sage: x = J.random_element()
1422 sage: while x.is_nilpotent() and not J.is_trivial():
1423 ....: x = J.random_element()
1424 sage: c = x.subalgebra_idempotent()
1429 if self
.parent().is_trivial():
1432 if self
.is_nilpotent():
1433 raise ValueError("this only works with non-nilpotent elements!")
1435 J
= self
.subalgebra_generated_by()
1438 # The image of the matrix of left-u^m-multiplication
1439 # will be minimal for some natural number s...
1441 minimal_dim
= J
.dimension()
1442 for i
in range(1, minimal_dim
):
1443 this_dim
= (u
**i
).operator().matrix().image().dimension()
1444 if this_dim
< minimal_dim
:
1445 minimal_dim
= this_dim
1448 # Now minimal_matrix should correspond to the smallest
1449 # non-zero subspace in Baes's (or really, Koecher's)
1452 # However, we need to restrict the matrix to work on the
1453 # subspace... or do we? Can't we just solve, knowing that
1454 # A(c) = u^(s+1) should have a solution in the big space,
1457 # Beware, solve_right() means that we're using COLUMN vectors.
1458 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1460 A
= u_next
.operator().matrix()
1461 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1463 # Now c is the idempotent we want, but it still lives in the subalgebra.
1464 return c
.superalgebra_element()
1469 Return my trace, the sum of my eigenvalues.
1471 In a trivial algebra, however you want to look at it, the trace is
1472 an empty sum for which we declare the result to be zero.
1476 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1483 sage: J = TrivialEJA()
1484 sage: J.zero().trace()
1488 sage: J = JordanSpinEJA(3)
1489 sage: x = sum(J.gens())
1495 sage: J = HadamardEJA(5)
1496 sage: J.one().trace()
1501 The trace of an element is a real number::
1503 sage: J = random_eja()
1504 sage: J.random_element().trace() in RLF
1507 The trace is linear::
1509 sage: J = random_eja()
1510 sage: x,y = J.random_elements(2)
1511 sage: alpha = J.base_ring().random_element()
1512 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1520 # Special case for the trivial algebra where
1521 # the trace is an empty sum.
1522 return P
.base_ring().zero()
1524 p
= P
._charpoly
_coefficients
()[r
-1]
1525 # The _charpoly_coeff function already adds the factor of
1526 # -1 to ensure that _charpoly_coeff(r-1) is really what
1527 # appears in front of t^{r-1} in the charpoly. However,
1528 # we want the negative of THAT for the trace.
1529 return -p(*self
.to_vector())
1532 def trace_inner_product(self
, other
):
1534 Return the trace inner product of myself and ``other``.
1538 sage: from mjo.eja.eja_algebra import random_eja
1542 The trace inner product is commutative, bilinear, and associative::
1544 sage: J = random_eja()
1545 sage: x,y,z = J.random_elements(3)
1547 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1550 sage: a = J.base_ring().random_element();
1551 sage: actual = (a*(x+z)).trace_inner_product(y)
1552 sage: expected = ( a*x.trace_inner_product(y) +
1553 ....: a*z.trace_inner_product(y) )
1554 sage: actual == expected
1556 sage: actual = x.trace_inner_product(a*(y+z))
1557 sage: expected = ( a*x.trace_inner_product(y) +
1558 ....: a*x.trace_inner_product(z) )
1559 sage: actual == expected
1562 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1566 if not other
in self
.parent():
1567 raise TypeError("'other' must live in the same algebra")
1569 return (self
*other
).trace()
1572 def trace_norm(self
):
1574 The norm of this element with respect to :meth:`trace_inner_product`.
1578 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1583 sage: J = HadamardEJA(2)
1584 sage: x = sum(J.gens())
1585 sage: x.trace_norm()
1590 sage: J = JordanSpinEJA(4)
1591 sage: x = sum(J.gens())
1592 sage: x.trace_norm()
1596 return self
.trace_inner_product(self
).sqrt()
1599 class CartesianProductEJAElement(FiniteDimensionalEJAElement
):
1602 Compute the determinant of this product-element using the
1603 determianants of its factors.
1605 This result Follows from the spectral decomposition of (say)
1606 the pair `(x,y)` in terms of the Jordan frame `\left\{ (c_1,
1607 0),(c_2, 0),...,(0,d_1),(0,d_2),... \right\}.
1609 from sage
.misc
.misc_c
import prod
1610 return prod( f
.det() for f
in self
.cartesian_factors() )
1612 def to_matrix(self
):
1613 # An override is necessary to call our custom _scale().
1614 B
= self
.parent().matrix_basis()
1615 W
= self
.parent().matrix_space()
1617 # Aaaaand linear combinations don't work in Cartesian
1618 # product spaces, even though they provide a method with
1619 # that name. This is hidden behind an "if" because the
1620 # _scale() function is slow.
1621 pairs
= zip(B
, self
.to_vector())
1622 return W
.sum( _scale(b
, alpha
) for (b
,alpha
) in pairs
)