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() # long time
488 sage: x = J.random_element() # long time
489 sage: (not x.operator().is_invertible()) or ( # long time
490 ....: x.operator().inverse()(J.one()) # long time
492 ....: x.inverse() ) # long time
495 Check that the fast (cached) and slow algorithms give the same
498 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
499 sage: x = J.random_element() # long time
500 sage: while not x.is_invertible(): # long time
501 ....: x = J.random_element() # long time
502 sage: slow = x.inverse() # long time
503 sage: _ = J._charpoly_coefficients() # long time
504 sage: fast = x.inverse() # long time
505 sage: slow == fast # long time
508 not_invertible_msg
= "element is not invertible"
509 if self
.parent()._charpoly
_coefficients
.is_in_cache():
510 # We can invert using our charpoly if it will be fast to
511 # compute. If the coefficients are cached, our rank had
513 if self
.det().is_zero():
514 raise ZeroDivisionError(not_invertible_msg
)
515 r
= self
.parent().rank()
516 a
= self
.characteristic_polynomial().coefficients(sparse
=False)
517 return (-1)**(r
+1)*sum(a
[i
+1]*self
**i
for i
in range(r
))/self
.det()
520 inv
= (~self
.quadratic_representation())(self
)
521 self
.is_invertible
.set_cache(True)
523 except ZeroDivisionError:
524 self
.is_invertible
.set_cache(False)
525 raise ZeroDivisionError(not_invertible_msg
)
529 def is_invertible(self
):
531 Return whether or not this element is invertible.
535 If computing my determinant will be fast, we do so and compare
536 with zero (Proposition II.2.4 in Faraut and
537 Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
538 reduces the problem to the invertibility of my quadratic
543 sage: from mjo.eja.eja_algebra import random_eja
547 The identity element is always invertible::
549 sage: J = random_eja()
550 sage: J.one().is_invertible()
553 The zero element is never invertible in a non-trivial algebra::
555 sage: J = random_eja()
556 sage: (not J.is_trivial()) and J.zero().is_invertible()
559 Test that the fast (cached) and slow algorithms give the same
562 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
563 sage: x = J.random_element() # long time
564 sage: slow = x.is_invertible() # long time
565 sage: _ = J._charpoly_coefficients() # long time
566 sage: fast = x.is_invertible() # long time
567 sage: slow == fast # long time
571 if self
.parent().is_trivial():
576 if self
.parent()._charpoly
_coefficients
.is_in_cache():
577 # The determinant will be quicker than inverting the
578 # quadratic representation, most likely.
579 return (not self
.det().is_zero())
581 # The easiest way to determine if I'm invertible is to try.
583 inv
= (~self
.quadratic_representation())(self
)
584 self
.inverse
.set_cache(inv
)
586 except ZeroDivisionError:
590 def is_primitive_idempotent(self
):
592 Return whether or not this element is a primitive (or minimal)
595 A primitive idempotent is a non-zero idempotent that is not
596 the sum of two other non-zero idempotents. Remark 2.7.15 in
597 Baes shows that this is what he refers to as a "minimal
600 An element of a Euclidean Jordan algebra is a minimal idempotent
601 if it :meth:`is_idempotent` and if its Peirce subalgebra
602 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
607 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
608 ....: RealSymmetricEJA,
614 This method is sloooooow.
618 The spectral decomposition of a non-regular element should always
619 contain at least one non-minimal idempotent::
621 sage: J = RealSymmetricEJA(3)
622 sage: x = sum(J.gens())
625 sage: [ c.is_primitive_idempotent()
626 ....: for (l,c) in x.spectral_decomposition() ]
629 On the other hand, the spectral decomposition of a regular
630 element should always be in terms of minimal idempotents::
632 sage: J = JordanSpinEJA(4)
633 sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
636 sage: [ c.is_primitive_idempotent()
637 ....: for (l,c) in x.spectral_decomposition() ]
642 The identity element is minimal only in an EJA of rank one::
644 sage: J = random_eja()
645 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
648 A non-idempotent cannot be a minimal idempotent::
650 sage: J = JordanSpinEJA(4)
651 sage: x = J.random_element()
652 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
655 Proposition 2.7.19 in Baes says that an element is a minimal
656 idempotent if and only if it's idempotent with trace equal to
659 sage: J = JordanSpinEJA(4)
660 sage: x = J.random_element()
661 sage: expected = (x.is_idempotent() and x.trace() == 1)
662 sage: actual = x.is_primitive_idempotent()
663 sage: actual == expected
666 Primitive idempotents must be non-zero::
668 sage: J = random_eja()
669 sage: J.zero().is_idempotent()
671 sage: J.zero().is_primitive_idempotent()
674 As a consequence of the fact that primitive idempotents must
675 be non-zero, there are no primitive idempotents in a trivial
676 Euclidean Jordan algebra::
678 sage: J = TrivialEJA()
679 sage: J.one().is_idempotent()
681 sage: J.one().is_primitive_idempotent()
685 if not self
.is_idempotent():
691 (_
,_
,J1
) = self
.parent().peirce_decomposition(self
)
692 return (J1
.dimension() == 1)
695 def is_nilpotent(self
):
697 Return whether or not some power of this element is zero.
701 We use Theorem 5 in Chapter III of Koecher, which says that
702 an element ``x`` is nilpotent if and only if ``x.operator()``
703 is nilpotent. And it is a basic fact of linear algebra that
704 an operator on an `n`-dimensional space is nilpotent if and
705 only if, when raised to the `n`th power, it equals the zero
706 operator (for example, see Axler Corollary 8.8).
710 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
715 sage: J = JordanSpinEJA(3)
716 sage: x = sum(J.gens())
717 sage: x.is_nilpotent()
722 The identity element is never nilpotent, except in a trivial EJA::
724 sage: J = random_eja()
725 sage: J.one().is_nilpotent() and not J.is_trivial()
728 The additive identity is always nilpotent::
730 sage: random_eja().zero().is_nilpotent()
735 zero_operator
= P
.zero().operator()
736 return self
.operator()**P
.dimension() == zero_operator
739 def is_regular(self
):
741 Return whether or not this is a regular element.
745 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
750 The identity element always has degree one, but any element
751 linearly-independent from it is regular::
753 sage: J = JordanSpinEJA(5)
754 sage: J.one().is_regular()
756 sage: b0, b1, b2, b3, b4 = J.gens()
759 sage: for x in J.gens():
760 ....: (J.one() + x).is_regular()
769 The zero element should never be regular, unless the parent
770 algebra has dimension less than or equal to one::
772 sage: J = random_eja()
773 sage: J.dimension() <= 1 or not J.zero().is_regular()
776 The unit element isn't regular unless the algebra happens to
777 consist of only its scalar multiples::
779 sage: J = random_eja()
780 sage: J.dimension() <= 1 or not J.one().is_regular()
784 return self
.degree() == self
.parent().rank()
789 Return the degree of this element, which is defined to be
790 the degree of its minimal polynomial.
798 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
803 sage: J = JordanSpinEJA(4)
804 sage: J.one().degree()
806 sage: b0,b1,b2,b3 = J.gens()
807 sage: (b0 - b1).degree()
810 In the spin factor algebra (of rank two), all elements that
811 aren't multiples of the identity are regular::
813 sage: J = JordanSpinEJA.random_instance()
814 sage: n = J.dimension()
815 sage: x = J.random_element()
816 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
821 The zero and unit elements are both of degree one in nontrivial
824 sage: J = random_eja()
825 sage: d = J.zero().degree()
826 sage: (J.is_trivial() and d == 0) or d == 1
828 sage: d = J.one().degree()
829 sage: (J.is_trivial() and d == 0) or d == 1
832 Our implementation agrees with the definition::
834 sage: x = random_eja().random_element()
835 sage: x.degree() == x.minimal_polynomial().degree()
839 n
= self
.parent().dimension()
842 # The minimal polynomial is an empty product, i.e. the
843 # constant polynomial "1" having degree zero.
846 # The minimal polynomial of zero in a nontrivial algebra
847 # is "t", and is of degree one.
850 # If this is a nonzero element of a nontrivial algebra, it
851 # has degree at least one. It follows that, in an algebra
852 # of dimension one, the degree must be actually one.
855 # BEWARE: The subalgebra_generated_by() method uses the result
856 # of this method to construct a basis for the subalgebra. That
857 # means, in particular, that we cannot implement this method
858 # as ``self.subalgebra_generated_by().dimension()``.
860 # Algorithm: keep appending (vector representations of) powers
861 # self as rows to a matrix and echelonizing it. When its rank
862 # stops increasing, we've reached a redundancy.
864 # Given the special cases above, we can assume that "self" is
865 # nonzero, the algebra is nontrivial, and that its dimension
867 M
= matrix([(self
.parent().one()).to_vector()])
870 # Specifying the row-reduction algorithm can e.g. help over
871 # AA because it avoids the RecursionError that gets thrown
872 # when we have to look too hard for a root.
874 # Beware: QQ supports an entirely different set of "algorithm"
875 # keywords than do AA and RR.
877 from sage
.rings
.all
import QQ
878 if self
.parent().base_ring() is not QQ
:
879 algo
= "scaled_partial_pivoting"
882 M
= matrix(M
.rows() + [(self
**d
).to_vector()])
885 if new_rank
== old_rank
:
894 def left_matrix(self
):
896 Our parent class defines ``left_matrix`` and ``matrix``
897 methods whose names are misleading. We don't want them.
899 raise NotImplementedError("use operator().matrix() instead")
904 def minimal_polynomial(self
):
906 Return the minimal polynomial of this element,
907 as a function of the variable `t`.
911 We restrict ourselves to the associative subalgebra
912 generated by this element, and then return the minimal
913 polynomial of this element's operator matrix (in that
914 subalgebra). This works by Baes Proposition 2.3.16.
918 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
919 ....: RealSymmetricEJA,
925 Keeping in mind that the polynomial ``1`` evaluates the identity
926 element (also the zero element) of the trivial algebra, it is clear
927 that the polynomial ``1`` is the minimal polynomial of the only
928 element in a trivial algebra::
930 sage: J = TrivialEJA()
931 sage: J.one().minimal_polynomial()
933 sage: J.zero().minimal_polynomial()
938 The minimal polynomial of the identity and zero elements are
939 always the same, except in trivial algebras where the minimal
940 polynomial of the unit/zero element is ``1``::
942 sage: J = random_eja()
943 sage: mu = J.one().minimal_polynomial()
944 sage: t = mu.parent().gen()
945 sage: mu + int(J.is_trivial())*(t-2)
947 sage: mu = J.zero().minimal_polynomial()
948 sage: t = mu.parent().gen()
949 sage: mu + int(J.is_trivial())*(t-1)
952 The degree of an element is (by one definition) the degree
953 of its minimal polynomial::
955 sage: x = random_eja().random_element()
956 sage: x.degree() == x.minimal_polynomial().degree()
959 The minimal polynomial and the characteristic polynomial coincide
960 and are known (see Alizadeh, Example 11.11) for all elements of
961 the spin factor algebra that aren't scalar multiples of the
962 identity. We require the dimension of the algebra to be at least
963 two here so that said elements actually exist::
965 sage: d_max = JordanSpinEJA._max_random_instance_dimension()
966 sage: n = ZZ.random_element(2, max(2,d_max))
967 sage: J = JordanSpinEJA(n)
968 sage: y = J.random_element()
969 sage: while y == y.coefficient(0)*J.one():
970 ....: y = J.random_element()
971 sage: y0 = y.to_vector()[0]
972 sage: y_bar = y.to_vector()[1:]
973 sage: actual = y.minimal_polynomial()
974 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
975 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
976 sage: bool(actual == expected)
979 The minimal polynomial should always kill its element::
981 sage: x = random_eja().random_element() # long time
982 sage: p = x.minimal_polynomial() # long time
983 sage: x.apply_univariate_polynomial(p) # long time
986 The minimal polynomial is invariant under a change of basis,
987 and in particular, a re-scaling of the basis::
989 sage: d_max = RealSymmetricEJA._max_random_instance_dimension()
990 sage: d = ZZ.random_element(1, d_max)
991 sage: n = RealSymmetricEJA._max_random_instance_size(d)
992 sage: J1 = RealSymmetricEJA(n)
993 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
994 sage: X = random_matrix(AA,n)
995 sage: X = X*X.transpose()
998 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
1003 # Pretty sure we know what the minimal polynomial of
1004 # the zero operator is going to be. This ensures
1005 # consistency of e.g. the polynomial variable returned
1006 # in the "normal" case without us having to think about it.
1007 return self
.operator().minimal_polynomial()
1009 # If we don't orthonormalize the subalgebra's basis, then the
1010 # first two monomials in the subalgebra will be self^0 and
1011 # self^1... assuming that self^1 is not a scalar multiple of
1012 # self^0 (the unit element). We special case these to avoid
1013 # having to solve a system to coerce self into the subalgebra.
1014 A
= self
.subalgebra_generated_by(orthonormalize
=False)
1016 if A
.dimension() == 1:
1017 # Does a solve to find the scalar multiple alpha such that
1018 # alpha*unit = self. We have to do this because the basis
1019 # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
1020 unit
= self
.parent().one()
1021 alpha
= self
.to_vector() / unit
.to_vector()
1022 return (unit
.operator()*alpha
).minimal_polynomial()
1024 # If the dimension of the subalgebra is >= 2, then we just
1025 # use the second basis element.
1026 return A
.monomial(1).operator().minimal_polynomial()
1030 def to_matrix(self
):
1032 Return an (often more natural) representation of this element as a
1035 Every finite-dimensional Euclidean Jordan Algebra is a direct
1036 sum of five simple algebras, four of which comprise Hermitian
1037 matrices. This method returns a "natural" matrix
1038 representation of this element as either a Hermitian matrix or
1043 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1045 ....: QuaternionHermitianEJA,
1046 ....: RealSymmetricEJA)
1050 sage: J = ComplexHermitianEJA(3)
1053 sage: J.one().to_matrix()
1064 sage: J = QuaternionHermitianEJA(2)
1067 sage: J.one().to_matrix()
1074 This also works in Cartesian product algebras::
1076 sage: J1 = HadamardEJA(1)
1077 sage: J2 = RealSymmetricEJA(2)
1078 sage: J = cartesian_product([J1,J2])
1079 sage: x = sum(J.gens())
1080 sage: x.to_matrix()[0]
1082 sage: x.to_matrix()[1]
1083 [ 1 0.7071067811865475?]
1084 [0.7071067811865475? 1]
1087 B
= self
.parent().matrix_basis()
1088 W
= self
.parent().matrix_space()
1090 # This is just a manual "from_vector()", but of course
1091 # matrix spaces aren't vector spaces in sage, so they
1092 # don't have a from_vector() method.
1093 return W
.linear_combination( zip(B
, self
.to_vector()) )
1099 The norm of this element with respect to :meth:`inner_product`.
1103 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1108 sage: J = HadamardEJA(2)
1109 sage: x = sum(J.gens())
1115 sage: J = JordanSpinEJA(4)
1116 sage: x = sum(J.gens())
1121 return self
.inner_product(self
).sqrt()
1126 Return the left-multiplication-by-this-element
1127 operator on the ambient algebra.
1131 sage: from mjo.eja.eja_algebra import random_eja
1135 sage: J = random_eja()
1136 sage: x,y = J.random_elements(2)
1137 sage: x.operator()(y) == x*y
1139 sage: y.operator()(x) == x*y
1144 left_mult_by_self
= lambda y
: self
*y
1145 L
= P
.module_morphism(function
=left_mult_by_self
, codomain
=P
)
1146 return FiniteDimensionalEJAOperator(P
, P
, L
.matrix() )
1149 def quadratic_representation(self
, other
=None):
1151 Return the quadratic representation of this element.
1155 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1160 The explicit form in the spin factor algebra is given by
1161 Alizadeh's Example 11.12::
1163 sage: x = JordanSpinEJA.random_instance().random_element()
1164 sage: x_vec = x.to_vector()
1165 sage: Q = matrix.identity(x.base_ring(), 0)
1166 sage: n = x_vec.degree()
1169 ....: x_bar = x_vec[1:]
1170 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1171 ....: B = 2*x0*x_bar.row()
1172 ....: C = 2*x0*x_bar.column()
1173 ....: D = matrix.identity(x.base_ring(), n-1)
1174 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1175 ....: D = D + 2*x_bar.tensor_product(x_bar)
1176 ....: Q = matrix.block(2,2,[A,B,C,D])
1177 sage: Q == x.quadratic_representation().matrix()
1180 Test all of the properties from Theorem 11.2 in Alizadeh::
1182 sage: J = random_eja()
1183 sage: x,y = J.random_elements(2)
1184 sage: Lx = x.operator()
1185 sage: Lxx = (x*x).operator()
1186 sage: Qx = x.quadratic_representation()
1187 sage: Qy = y.quadratic_representation()
1188 sage: Qxy = x.quadratic_representation(y)
1189 sage: Qex = J.one().quadratic_representation(x)
1190 sage: n = ZZ.random_element(10)
1191 sage: Qxn = (x^n).quadratic_representation()
1195 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1198 Property 2 (multiply on the right for :trac:`28272`):
1200 sage: alpha = J.base_ring().random_element()
1201 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1206 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1209 sage: not x.is_invertible() or (
1212 ....: x.inverse().quadratic_representation() )
1215 sage: Qxy(J.one()) == x*y
1220 sage: not x.is_invertible() or (
1221 ....: x.quadratic_representation(x.inverse())*Qx
1222 ....: == Qx*x.quadratic_representation(x.inverse()) )
1225 sage: not x.is_invertible() or (
1226 ....: x.quadratic_representation(x.inverse())*Qx
1228 ....: 2*Lx*Qex - Qx )
1231 sage: 2*Lx*Qex - Qx == Lxx
1236 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1246 sage: not x.is_invertible() or (
1247 ....: Qx*x.inverse().operator() == Lx )
1252 sage: not x.operator_commutes_with(y) or (
1253 ....: Qx(y)^n == Qxn(y^n) )
1259 elif not other
in self
.parent():
1260 raise TypeError("'other' must live in the same algebra")
1263 M
= other
.operator()
1264 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1268 def spectral_decomposition(self
):
1270 Return the unique spectral decomposition of this element.
1274 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1275 element's left-multiplication-by operator to the subalgebra it
1276 generates. We then compute the spectral decomposition of that
1277 operator, and the spectral projectors we get back must be the
1278 left-multiplication-by operators for the idempotents we
1279 seek. Thus applying them to the identity element gives us those
1282 Since the eigenvalues are required to be distinct, we take
1283 the spectral decomposition of the zero element to be zero
1284 times the identity element of the algebra (which is idempotent,
1289 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1293 The spectral decomposition of the identity is ``1`` times itself,
1294 and the spectral decomposition of zero is ``0`` times the identity::
1296 sage: J = RealSymmetricEJA(3)
1299 sage: J.one().spectral_decomposition()
1301 sage: J.zero().spectral_decomposition()
1306 sage: J = RealSymmetricEJA(4)
1307 sage: x = sum(J.gens())
1308 sage: sd = x.spectral_decomposition()
1313 sage: c0.inner_product(c1) == 0
1315 sage: c0.is_idempotent()
1317 sage: c1.is_idempotent()
1319 sage: c0 + c1 == J.one()
1321 sage: l0*c0 + l1*c1 == x
1324 The spectral decomposition should work in subalgebras, too::
1326 sage: J = RealSymmetricEJA(4)
1327 sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
1328 sage: A = 2*b5 - 2*b8
1329 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1330 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1331 sage: (f0, f1, f2) = J1.gens()
1332 sage: f0.spectral_decomposition()
1333 [(0, 1.000000000000000?*c2), (1, 1.000000000000000?*c0)]
1336 A
= self
.subalgebra_generated_by(orthonormalize
=True)
1338 for (evalue
, proj
) in A(self
).operator().spectral_decomposition():
1339 result
.append( (evalue
, proj(A
.one()).superalgebra_element()) )
1342 def subalgebra_generated_by(self
, **kwargs
):
1344 Return the associative subalgebra of the parent EJA generated
1347 Since our parent algebra is unital, we want "subalgebra" to mean
1348 "unital subalgebra" as well; thus the subalgebra that an element
1349 generates will itself be a Euclidean Jordan algebra after
1350 restricting the algebra operations appropriately. This is the
1351 subalgebra that Faraut and Korányi work with in section II.2, for
1356 sage: from mjo.eja.eja_algebra import (random_eja,
1358 ....: RealSymmetricEJA)
1362 We can create subalgebras of Cartesian product EJAs that are not
1363 themselves Cartesian product EJAs (they're just "regular" EJAs)::
1365 sage: J1 = HadamardEJA(3)
1366 sage: J2 = RealSymmetricEJA(2)
1367 sage: J = cartesian_product([J1,J2])
1368 sage: J.one().subalgebra_generated_by()
1369 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
1373 This subalgebra, being composed of only powers, is associative::
1375 sage: x0 = random_eja().random_element()
1376 sage: A = x0.subalgebra_generated_by()
1377 sage: x,y,z = A.random_elements(3)
1378 sage: (x*y)*z == x*(y*z)
1381 Squaring in the subalgebra should work the same as in
1384 sage: x = random_eja().random_element()
1385 sage: A = x.subalgebra_generated_by()
1386 sage: A(x^2) == A(x)*A(x)
1389 By definition, the subalgebra generated by the zero element is
1390 the one-dimensional algebra generated by the identity
1391 element... unless the original algebra was trivial, in which
1392 case the subalgebra is trivial too::
1394 sage: A = random_eja().zero().subalgebra_generated_by()
1395 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1399 powers
= tuple( self
**k
for k
in range(self
.degree()) )
1400 A
= self
.parent().subalgebra(powers
,
1405 A
.one
.set_cache(A(self
.parent().one()))
1409 def subalgebra_idempotent(self
):
1411 Find an idempotent in the associative subalgebra I generate
1412 using Proposition 2.3.5 in Baes.
1416 sage: from mjo.eja.eja_algebra import random_eja
1420 Ensure that we can find an idempotent in a non-trivial algebra
1421 where there are non-nilpotent elements, or that we get the dumb
1422 solution in the trivial algebra::
1424 sage: J = random_eja()
1425 sage: x = J.random_element()
1426 sage: while x.is_nilpotent() and not J.is_trivial():
1427 ....: x = J.random_element()
1428 sage: c = x.subalgebra_idempotent()
1433 if self
.parent().is_trivial():
1436 if self
.is_nilpotent():
1437 raise ValueError("this only works with non-nilpotent elements!")
1439 J
= self
.subalgebra_generated_by()
1442 # The image of the matrix of left-u^m-multiplication
1443 # will be minimal for some natural number s...
1445 minimal_dim
= J
.dimension()
1446 for i
in range(1, minimal_dim
):
1447 this_dim
= (u
**i
).operator().matrix().image().dimension()
1448 if this_dim
< minimal_dim
:
1449 minimal_dim
= this_dim
1452 # Now minimal_matrix should correspond to the smallest
1453 # non-zero subspace in Baes's (or really, Koecher's)
1456 # However, we need to restrict the matrix to work on the
1457 # subspace... or do we? Can't we just solve, knowing that
1458 # A(c) = u^(s+1) should have a solution in the big space,
1461 # Beware, solve_right() means that we're using COLUMN vectors.
1462 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1464 A
= u_next
.operator().matrix()
1465 c
= J
.from_vector(A
.solve_right(u_next
.to_vector()))
1467 # Now c is the idempotent we want, but it still lives in the subalgebra.
1468 return c
.superalgebra_element()
1473 Return my trace, the sum of my eigenvalues.
1475 In a trivial algebra, however you want to look at it, the trace is
1476 an empty sum for which we declare the result to be zero.
1480 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1487 sage: J = TrivialEJA()
1488 sage: J.zero().trace()
1492 sage: J = JordanSpinEJA(3)
1493 sage: x = sum(J.gens())
1499 sage: J = HadamardEJA(5)
1500 sage: J.one().trace()
1505 The trace of an element is a real number::
1507 sage: J = random_eja()
1508 sage: J.random_element().trace() in RLF
1511 The trace is linear::
1513 sage: J = random_eja()
1514 sage: x,y = J.random_elements(2)
1515 sage: alpha = J.base_ring().random_element()
1516 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1519 The trace of a square is nonnegative::
1521 sage: x = random_eja().random_element()
1522 sage: (x*x).trace() >= 0
1530 # Special case for the trivial algebra where
1531 # the trace is an empty sum.
1532 return P
.base_ring().zero()
1534 p
= P
._charpoly
_coefficients
()[r
-1]
1535 # The _charpoly_coeff function already adds the factor of
1536 # -1 to ensure that _charpoly_coeff(r-1) is really what
1537 # appears in front of t^{r-1} in the charpoly. However,
1538 # we want the negative of THAT for the trace.
1539 return -p(*self
.to_vector())
1542 def trace_inner_product(self
, other
):
1544 Return the trace inner product of myself and ``other``.
1548 sage: from mjo.eja.eja_algebra import random_eja
1552 The trace inner product is commutative, bilinear, and associative::
1554 sage: J = random_eja()
1555 sage: x,y,z = J.random_elements(3)
1557 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1560 sage: a = J.base_ring().random_element();
1561 sage: actual = (a*(x+z)).trace_inner_product(y)
1562 sage: expected = ( a*x.trace_inner_product(y) +
1563 ....: a*z.trace_inner_product(y) )
1564 sage: actual == expected
1566 sage: actual = x.trace_inner_product(a*(y+z))
1567 sage: expected = ( a*x.trace_inner_product(y) +
1568 ....: a*x.trace_inner_product(z) )
1569 sage: actual == expected
1572 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1576 if not other
in self
.parent():
1577 raise TypeError("'other' must live in the same algebra")
1579 return (self
*other
).trace()
1582 def trace_norm(self
):
1584 The norm of this element with respect to :meth:`trace_inner_product`.
1588 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1593 sage: J = HadamardEJA(2)
1594 sage: x = sum(J.gens())
1595 sage: x.trace_norm()
1600 sage: J = JordanSpinEJA(4)
1601 sage: x = sum(J.gens())
1602 sage: x.trace_norm()
1606 return self
.trace_inner_product(self
).sqrt()
1609 class CartesianProductEJAElement(FiniteDimensionalEJAElement
):
1612 Compute the determinant of this product-element using the
1613 determianants of its factors.
1615 This result Follows from the spectral decomposition of (say)
1616 the pair `(x,y)` in terms of the Jordan frame `\left\{ (c_1,
1617 0),(c_2, 0),...,(0,d_1),(0,d_2),... \right\}.
1619 from sage
.misc
.misc_c
import prod
1620 return prod( f
.det() for f
in self
.cartesian_factors() )
1622 def to_matrix(self
):
1623 # An override is necessary to call our custom _scale().
1624 B
= self
.parent().matrix_basis()
1625 W
= self
.parent().matrix_space()
1627 # Aaaaand linear combinations don't work in Cartesian
1628 # product spaces, even though they provide a method with
1629 # that name. This is hidden behind an "if" because the
1630 # _scale() function is slow.
1631 pairs
= zip(B
, self
.to_vector())
1632 return W
.sum( _scale(b
, alpha
) for (b
,alpha
) in pairs
)