2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
9 from sage
.structure
.element
import is_Matrix
10 from sage
.structure
.category_object
import normalize_names
12 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra
import FiniteDimensionalAlgebra
13 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra_element
import FiniteDimensionalAlgebraElement
15 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra
):
17 def __classcall_private__(cls
,
21 assume_associative
=False,
26 mult_table
= [b
.base_extend(field
) for b
in mult_table
]
29 if not (is_Matrix(b
) and b
.dimensions() == (n
, n
)):
30 raise ValueError("input is not a multiplication table")
31 mult_table
= tuple(mult_table
)
33 cat
= MagmaticAlgebras(field
).FiniteDimensional().WithBasis()
34 cat
.or_subcategory(category
)
35 if assume_associative
:
36 cat
= cat
.Associative()
38 names
= normalize_names(n
, names
)
40 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, cls
)
41 return fda
.__classcall
__(cls
,
44 assume_associative
=assume_associative
,
48 natural_basis
=natural_basis
)
55 assume_associative
=False,
62 By definition, Jordan multiplication commutes::
64 sage: set_random_seed()
65 sage: J = random_eja()
66 sage: x = J.random_element()
67 sage: y = J.random_element()
73 self
._natural
_basis
= natural_basis
74 self
._multiplication
_table
= mult_table
75 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
84 Return a string representation of ``self``.
86 fmt
= "Euclidean Jordan algebra of degree {} over {}"
87 return fmt
.format(self
.degree(), self
.base_ring())
90 def inner_product(self
, x
, y
):
92 The inner product associated with this Euclidean Jordan algebra.
94 Defaults to the trace inner product, but can be overridden by
95 subclasses if they are sure that the necessary properties are
100 The inner product must satisfy its axiom for this algebra to truly
101 be a Euclidean Jordan Algebra::
103 sage: set_random_seed()
104 sage: J = random_eja()
105 sage: x = J.random_element()
106 sage: y = J.random_element()
107 sage: z = J.random_element()
108 sage: (x*y).inner_product(z) == y.inner_product(x*z)
112 if (not x
in self
) or (not y
in self
):
113 raise TypeError("arguments must live in this algebra")
114 return x
.trace_inner_product(y
)
117 def natural_basis(self
):
119 Return a more-natural representation of this algebra's basis.
121 Every finite-dimensional Euclidean Jordan Algebra is a direct
122 sum of five simple algebras, four of which comprise Hermitian
123 matrices. This method returns the original "natural" basis
124 for our underlying vector space. (Typically, the natural basis
125 is used to construct the multiplication table in the first place.)
127 Note that this will always return a matrix. The standard basis
128 in `R^n` will be returned as `n`-by-`1` column matrices.
132 sage: J = RealSymmetricEJA(2)
135 sage: J.natural_basis()
143 sage: J = JordanSpinEJA(2)
146 sage: J.natural_basis()
153 if self
._natural
_basis
is None:
154 return tuple( b
.vector().column() for b
in self
.basis() )
156 return self
._natural
_basis
161 Return the rank of this EJA.
163 if self
._rank
is None:
164 raise ValueError("no rank specified at genesis")
169 class Element(FiniteDimensionalAlgebraElement
):
171 An element of a Euclidean Jordan algebra.
174 def __init__(self
, A
, elt
=None):
178 The identity in `S^n` is converted to the identity in the EJA::
180 sage: J = RealSymmetricEJA(3)
181 sage: I = identity_matrix(QQ,3)
182 sage: J(I) == J.one()
185 This skew-symmetric matrix can't be represented in the EJA::
187 sage: J = RealSymmetricEJA(3)
188 sage: A = matrix(QQ,3, lambda i,j: i-j)
190 Traceback (most recent call last):
192 ArithmeticError: vector is not in free module
195 # Goal: if we're given a matrix, and if it lives in our
196 # parent algebra's "natural ambient space," convert it
197 # into an algebra element.
199 # The catch is, we make a recursive call after converting
200 # the given matrix into a vector that lives in the algebra.
201 # This we need to try the parent class initializer first,
202 # to avoid recursing forever if we're given something that
203 # already fits into the algebra, but also happens to live
204 # in the parent's "natural ambient space" (this happens with
207 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, elt
)
209 natural_basis
= A
.natural_basis()
210 if elt
in natural_basis
[0].matrix_space():
211 # Thanks for nothing! Matrix spaces aren't vector
212 # spaces in Sage, so we have to figure out its
213 # natural-basis coordinates ourselves.
214 V
= VectorSpace(elt
.base_ring(), elt
.nrows()**2)
215 W
= V
.span( _mat2vec(s
) for s
in natural_basis
)
216 coords
= W
.coordinates(_mat2vec(elt
))
217 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, coords
)
219 def __pow__(self
, n
):
221 Return ``self`` raised to the power ``n``.
223 Jordan algebras are always power-associative; see for
224 example Faraut and Koranyi, Proposition II.1.2 (ii).
228 We have to override this because our superclass uses row vectors
229 instead of column vectors! We, on the other hand, assume column
234 sage: set_random_seed()
235 sage: x = random_eja().random_element()
236 sage: x.operator_matrix()*x.vector() == (x^2).vector()
239 A few examples of power-associativity::
241 sage: set_random_seed()
242 sage: x = random_eja().random_element()
243 sage: x*(x*x)*(x*x) == x^5
245 sage: (x*x)*(x*x*x) == x^5
248 We also know that powers operator-commute (Koecher, Chapter
251 sage: set_random_seed()
252 sage: x = random_eja().random_element()
253 sage: m = ZZ.random_element(0,10)
254 sage: n = ZZ.random_element(0,10)
255 sage: Lxm = (x^m).operator_matrix()
256 sage: Lxn = (x^n).operator_matrix()
257 sage: Lxm*Lxn == Lxn*Lxm
267 return A( (self
.operator_matrix()**(n
-1))*self
.vector() )
270 def characteristic_polynomial(self
):
272 Return my characteristic polynomial (if I'm a regular
275 Eventually this should be implemented in terms of the parent
276 algebra's characteristic polynomial that works for ALL
279 if self
.is_regular():
280 return self
.minimal_polynomial()
282 raise NotImplementedError('irregular element')
285 def inner_product(self
, other
):
287 Return the parent algebra's inner product of myself and ``other``.
291 The inner product in the Jordan spin algebra is the usual
292 inner product on `R^n` (this example only works because the
293 basis for the Jordan algebra is the standard basis in `R^n`)::
295 sage: J = JordanSpinEJA(3)
296 sage: x = vector(QQ,[1,2,3])
297 sage: y = vector(QQ,[4,5,6])
298 sage: x.inner_product(y)
300 sage: J(x).inner_product(J(y))
303 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
304 multiplication is the usual matrix multiplication in `S^n`,
305 so the inner product of the identity matrix with itself
308 sage: J = RealSymmetricEJA(3)
309 sage: J.one().inner_product(J.one())
312 Likewise, the inner product on `C^n` is `<X,Y> =
313 Re(trace(X*Y))`, where we must necessarily take the real
314 part because the product of Hermitian matrices may not be
317 sage: J = ComplexHermitianEJA(3)
318 sage: J.one().inner_product(J.one())
321 Ditto for the quaternions::
323 sage: J = QuaternionHermitianEJA(3)
324 sage: J.one().inner_product(J.one())
329 Ensure that we can always compute an inner product, and that
330 it gives us back a real number::
332 sage: set_random_seed()
333 sage: J = random_eja()
334 sage: x = J.random_element()
335 sage: y = J.random_element()
336 sage: x.inner_product(y) in RR
342 raise TypeError("'other' must live in the same algebra")
344 return P
.inner_product(self
, other
)
347 def operator_commutes_with(self
, other
):
349 Return whether or not this element operator-commutes
354 The definition of a Jordan algebra says that any element
355 operator-commutes with its square::
357 sage: set_random_seed()
358 sage: x = random_eja().random_element()
359 sage: x.operator_commutes_with(x^2)
364 Test Lemma 1 from Chapter III of Koecher::
366 sage: set_random_seed()
367 sage: J = random_eja()
368 sage: u = J.random_element()
369 sage: v = J.random_element()
370 sage: lhs = u.operator_commutes_with(u*v)
371 sage: rhs = v.operator_commutes_with(u^2)
376 if not other
in self
.parent():
377 raise TypeError("'other' must live in the same algebra")
379 A
= self
.operator_matrix()
380 B
= other
.operator_matrix()
386 Return my determinant, the product of my eigenvalues.
390 sage: J = JordanSpinEJA(2)
391 sage: e0,e1 = J.gens()
395 sage: J = JordanSpinEJA(3)
396 sage: e0,e1,e2 = J.gens()
397 sage: x = e0 + e1 + e2
402 cs
= self
.characteristic_polynomial().coefficients(sparse
=False)
405 return cs
[0] * (-1)**r
407 raise ValueError('charpoly had no coefficients')
412 Return the Jordan-multiplicative inverse of this element.
414 We can't use the superclass method because it relies on the
415 algebra being associative.
419 The inverse in the spin factor algebra is given in Alizadeh's
422 sage: set_random_seed()
423 sage: n = ZZ.random_element(1,10)
424 sage: J = JordanSpinEJA(n)
425 sage: x = J.random_element()
426 sage: while x.is_zero():
427 ....: x = J.random_element()
428 sage: x_vec = x.vector()
430 sage: x_bar = x_vec[1:]
431 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
432 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
433 sage: x_inverse = coeff*inv_vec
434 sage: x.inverse() == J(x_inverse)
439 The identity element is its own inverse::
441 sage: set_random_seed()
442 sage: J = random_eja()
443 sage: J.one().inverse() == J.one()
446 If an element has an inverse, it acts like one. TODO: this
447 can be a lot less ugly once ``is_invertible`` doesn't crash
448 on irregular elements::
450 sage: set_random_seed()
451 sage: J = random_eja()
452 sage: x = J.random_element()
454 ....: x.inverse()*x == J.one()
460 if self
.parent().is_associative():
461 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
464 # TODO: we can do better once the call to is_invertible()
465 # doesn't crash on irregular elements.
466 #if not self.is_invertible():
467 # raise ValueError('element is not invertible')
469 # We do this a little different than the usual recursive
470 # call to a finite-dimensional algebra element, because we
471 # wind up with an inverse that lives in the subalgebra and
472 # we need information about the parent to convert it back.
473 V
= self
.span_of_powers()
474 assoc_subalg
= self
.subalgebra_generated_by()
475 # Mis-design warning: the basis used for span_of_powers()
476 # and subalgebra_generated_by() must be the same, and in
478 elt
= assoc_subalg(V
.coordinates(self
.vector()))
480 # This will be in the subalgebra's coordinates...
481 fda_elt
= FiniteDimensionalAlgebraElement(assoc_subalg
, elt
)
482 subalg_inverse
= fda_elt
.inverse()
484 # So we have to convert back...
485 basis
= [ self
.parent(v
) for v
in V
.basis() ]
486 pairs
= zip(subalg_inverse
.vector(), basis
)
487 return self
.parent().linear_combination(pairs
)
490 def is_invertible(self
):
492 Return whether or not this element is invertible.
494 We can't use the superclass method because it relies on
495 the algebra being associative.
499 The usual way to do this is to check if the determinant is
500 zero, but we need the characteristic polynomial for the
501 determinant. The minimal polynomial is a lot easier to get,
502 so we use Corollary 2 in Chapter V of Koecher to check
503 whether or not the paren't algebra's zero element is a root
504 of this element's minimal polynomial.
508 The identity element is always invertible::
510 sage: set_random_seed()
511 sage: J = random_eja()
512 sage: J.one().is_invertible()
515 The zero element is never invertible::
517 sage: set_random_seed()
518 sage: J = random_eja()
519 sage: J.zero().is_invertible()
523 zero
= self
.parent().zero()
524 p
= self
.minimal_polynomial()
525 return not (p(zero
) == zero
)
528 def is_nilpotent(self
):
530 Return whether or not some power of this element is zero.
532 The superclass method won't work unless we're in an
533 associative algebra, and we aren't. However, we generate
534 an assocoative subalgebra and we're nilpotent there if and
535 only if we're nilpotent here (probably).
539 The identity element is never nilpotent::
541 sage: set_random_seed()
542 sage: random_eja().one().is_nilpotent()
545 The additive identity is always nilpotent::
547 sage: set_random_seed()
548 sage: random_eja().zero().is_nilpotent()
552 # The element we're going to call "is_nilpotent()" on.
553 # Either myself, interpreted as an element of a finite-
554 # dimensional algebra, or an element of an associative
558 if self
.parent().is_associative():
559 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
561 V
= self
.span_of_powers()
562 assoc_subalg
= self
.subalgebra_generated_by()
563 # Mis-design warning: the basis used for span_of_powers()
564 # and subalgebra_generated_by() must be the same, and in
566 elt
= assoc_subalg(V
.coordinates(self
.vector()))
568 # Recursive call, but should work since elt lives in an
569 # associative algebra.
570 return elt
.is_nilpotent()
573 def is_regular(self
):
575 Return whether or not this is a regular element.
579 The identity element always has degree one, but any element
580 linearly-independent from it is regular::
582 sage: J = JordanSpinEJA(5)
583 sage: J.one().is_regular()
585 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
586 sage: for x in J.gens():
587 ....: (J.one() + x).is_regular()
595 return self
.degree() == self
.parent().rank()
600 Compute the degree of this element the straightforward way
601 according to the definition; by appending powers to a list
602 and figuring out its dimension (that is, whether or not
603 they're linearly dependent).
607 sage: J = JordanSpinEJA(4)
608 sage: J.one().degree()
610 sage: e0,e1,e2,e3 = J.gens()
611 sage: (e0 - e1).degree()
614 In the spin factor algebra (of rank two), all elements that
615 aren't multiples of the identity are regular::
617 sage: set_random_seed()
618 sage: n = ZZ.random_element(1,10)
619 sage: J = JordanSpinEJA(n)
620 sage: x = J.random_element()
621 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
625 return self
.span_of_powers().dimension()
628 def minimal_polynomial(self
):
632 We restrict ourselves to the associative subalgebra
633 generated by this element, and then return the minimal
634 polynomial of this element's operator matrix (in that
635 subalgebra). This works by Baes Proposition 2.3.16.
639 sage: set_random_seed()
640 sage: x = random_eja().random_element()
641 sage: x.degree() == x.minimal_polynomial().degree()
646 sage: set_random_seed()
647 sage: x = random_eja().random_element()
648 sage: x.degree() == x.minimal_polynomial().degree()
651 The minimal polynomial and the characteristic polynomial coincide
652 and are known (see Alizadeh, Example 11.11) for all elements of
653 the spin factor algebra that aren't scalar multiples of the
656 sage: set_random_seed()
657 sage: n = ZZ.random_element(2,10)
658 sage: J = JordanSpinEJA(n)
659 sage: y = J.random_element()
660 sage: while y == y.coefficient(0)*J.one():
661 ....: y = J.random_element()
662 sage: y0 = y.vector()[0]
663 sage: y_bar = y.vector()[1:]
664 sage: actual = y.minimal_polynomial()
665 sage: x = SR.symbol('x', domain='real')
666 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
667 sage: bool(actual == expected)
671 V
= self
.span_of_powers()
672 assoc_subalg
= self
.subalgebra_generated_by()
673 # Mis-design warning: the basis used for span_of_powers()
674 # and subalgebra_generated_by() must be the same, and in
676 elt
= assoc_subalg(V
.coordinates(self
.vector()))
677 return elt
.operator_matrix().minimal_polynomial()
680 def natural_representation(self
):
682 Return a more-natural representation of this element.
684 Every finite-dimensional Euclidean Jordan Algebra is a
685 direct sum of five simple algebras, four of which comprise
686 Hermitian matrices. This method returns the original
687 "natural" representation of this element as a Hermitian
688 matrix, if it has one. If not, you get the usual representation.
692 sage: J = ComplexHermitianEJA(3)
695 sage: J.one().natural_representation()
705 sage: J = QuaternionHermitianEJA(3)
708 sage: J.one().natural_representation()
709 [1 0 0 0 0 0 0 0 0 0 0 0]
710 [0 1 0 0 0 0 0 0 0 0 0 0]
711 [0 0 1 0 0 0 0 0 0 0 0 0]
712 [0 0 0 1 0 0 0 0 0 0 0 0]
713 [0 0 0 0 1 0 0 0 0 0 0 0]
714 [0 0 0 0 0 1 0 0 0 0 0 0]
715 [0 0 0 0 0 0 1 0 0 0 0 0]
716 [0 0 0 0 0 0 0 1 0 0 0 0]
717 [0 0 0 0 0 0 0 0 1 0 0 0]
718 [0 0 0 0 0 0 0 0 0 1 0 0]
719 [0 0 0 0 0 0 0 0 0 0 1 0]
720 [0 0 0 0 0 0 0 0 0 0 0 1]
723 B
= self
.parent().natural_basis()
724 W
= B
[0].matrix_space()
725 return W
.linear_combination(zip(self
.vector(), B
))
728 def operator_matrix(self
):
730 Return the matrix that represents left- (or right-)
731 multiplication by this element in the parent algebra.
733 We have to override this because the superclass method
734 returns a matrix that acts on row vectors (that is, on
739 Test the first polarization identity from my notes, Koecher Chapter
740 III, or from Baes (2.3)::
742 sage: set_random_seed()
743 sage: J = random_eja()
744 sage: x = J.random_element()
745 sage: y = J.random_element()
746 sage: Lx = x.operator_matrix()
747 sage: Ly = y.operator_matrix()
748 sage: Lxx = (x*x).operator_matrix()
749 sage: Lxy = (x*y).operator_matrix()
750 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
753 Test the second polarization identity from my notes or from
756 sage: set_random_seed()
757 sage: J = random_eja()
758 sage: x = J.random_element()
759 sage: y = J.random_element()
760 sage: z = J.random_element()
761 sage: Lx = x.operator_matrix()
762 sage: Ly = y.operator_matrix()
763 sage: Lz = z.operator_matrix()
764 sage: Lzy = (z*y).operator_matrix()
765 sage: Lxy = (x*y).operator_matrix()
766 sage: Lxz = (x*z).operator_matrix()
767 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
770 Test the third polarization identity from my notes or from
773 sage: set_random_seed()
774 sage: J = random_eja()
775 sage: u = J.random_element()
776 sage: y = J.random_element()
777 sage: z = J.random_element()
778 sage: Lu = u.operator_matrix()
779 sage: Ly = y.operator_matrix()
780 sage: Lz = z.operator_matrix()
781 sage: Lzy = (z*y).operator_matrix()
782 sage: Luy = (u*y).operator_matrix()
783 sage: Luz = (u*z).operator_matrix()
784 sage: Luyz = (u*(y*z)).operator_matrix()
785 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
786 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
787 sage: bool(lhs == rhs)
791 fda_elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
792 return fda_elt
.matrix().transpose()
795 def quadratic_representation(self
, other
=None):
797 Return the quadratic representation of this element.
801 The explicit form in the spin factor algebra is given by
802 Alizadeh's Example 11.12::
804 sage: set_random_seed()
805 sage: n = ZZ.random_element(1,10)
806 sage: J = JordanSpinEJA(n)
807 sage: x = J.random_element()
808 sage: x_vec = x.vector()
810 sage: x_bar = x_vec[1:]
811 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
812 sage: B = 2*x0*x_bar.row()
813 sage: C = 2*x0*x_bar.column()
814 sage: D = identity_matrix(QQ, n-1)
815 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
816 sage: D = D + 2*x_bar.tensor_product(x_bar)
817 sage: Q = block_matrix(2,2,[A,B,C,D])
818 sage: Q == x.quadratic_representation()
821 Test all of the properties from Theorem 11.2 in Alizadeh::
823 sage: set_random_seed()
824 sage: J = random_eja()
825 sage: x = J.random_element()
826 sage: y = J.random_element()
830 sage: actual = x.quadratic_representation(y)
831 sage: expected = ( (x+y).quadratic_representation()
832 ....: -x.quadratic_representation()
833 ....: -y.quadratic_representation() ) / 2
834 sage: actual == expected
839 sage: alpha = QQ.random_element()
840 sage: actual = (alpha*x).quadratic_representation()
841 sage: expected = (alpha^2)*x.quadratic_representation()
842 sage: actual == expected
847 sage: Qy = y.quadratic_representation()
848 sage: actual = J(Qy*x.vector()).quadratic_representation()
849 sage: expected = Qy*x.quadratic_representation()*Qy
850 sage: actual == expected
855 sage: k = ZZ.random_element(1,10)
856 sage: actual = (x^k).quadratic_representation()
857 sage: expected = (x.quadratic_representation())^k
858 sage: actual == expected
864 elif not other
in self
.parent():
865 raise TypeError("'other' must live in the same algebra")
867 L
= self
.operator_matrix()
868 M
= other
.operator_matrix()
869 return ( L
*M
+ M
*L
- (self
*other
).operator_matrix() )
872 def span_of_powers(self
):
874 Return the vector space spanned by successive powers of
877 # The dimension of the subalgebra can't be greater than
878 # the big algebra, so just put everything into a list
879 # and let span() get rid of the excess.
881 # We do the extra ambient_vector_space() in case we're messing
882 # with polynomials and the direct parent is a module.
883 V
= self
.vector().parent().ambient_vector_space()
884 return V
.span( (self
**d
).vector() for d
in xrange(V
.dimension()) )
887 def subalgebra_generated_by(self
):
889 Return the associative subalgebra of the parent EJA generated
894 sage: set_random_seed()
895 sage: x = random_eja().random_element()
896 sage: x.subalgebra_generated_by().is_associative()
899 Squaring in the subalgebra should be the same thing as
900 squaring in the superalgebra::
902 sage: set_random_seed()
903 sage: x = random_eja().random_element()
904 sage: u = x.subalgebra_generated_by().random_element()
905 sage: u.operator_matrix()*u.vector() == (u**2).vector()
909 # First get the subspace spanned by the powers of myself...
910 V
= self
.span_of_powers()
913 # Now figure out the entries of the right-multiplication
914 # matrix for the successive basis elements b0, b1,... of
917 for b_right
in V
.basis():
918 eja_b_right
= self
.parent()(b_right
)
920 # The first row of the right-multiplication matrix by
921 # b1 is what we get if we apply that matrix to b1. The
922 # second row of the right multiplication matrix by b1
923 # is what we get when we apply that matrix to b2...
925 # IMPORTANT: this assumes that all vectors are COLUMN
926 # vectors, unlike our superclass (which uses row vectors).
927 for b_left
in V
.basis():
928 eja_b_left
= self
.parent()(b_left
)
929 # Multiply in the original EJA, but then get the
930 # coordinates from the subalgebra in terms of its
932 this_row
= V
.coordinates((eja_b_left
*eja_b_right
).vector())
933 b_right_rows
.append(this_row
)
934 b_right_matrix
= matrix(F
, b_right_rows
)
935 mats
.append(b_right_matrix
)
937 # It's an algebra of polynomials in one element, and EJAs
938 # are power-associative.
940 # TODO: choose generator names intelligently.
941 return FiniteDimensionalEuclideanJordanAlgebra(F
, mats
, assume_associative
=True, names
='f')
944 def subalgebra_idempotent(self
):
946 Find an idempotent in the associative subalgebra I generate
947 using Proposition 2.3.5 in Baes.
951 sage: set_random_seed()
952 sage: J = RealCartesianProductEJA(5)
953 sage: c = J.random_element().subalgebra_idempotent()
956 sage: J = JordanSpinEJA(5)
957 sage: c = J.random_element().subalgebra_idempotent()
962 if self
.is_nilpotent():
963 raise ValueError("this only works with non-nilpotent elements!")
965 V
= self
.span_of_powers()
966 J
= self
.subalgebra_generated_by()
967 # Mis-design warning: the basis used for span_of_powers()
968 # and subalgebra_generated_by() must be the same, and in
970 u
= J(V
.coordinates(self
.vector()))
972 # The image of the matrix of left-u^m-multiplication
973 # will be minimal for some natural number s...
975 minimal_dim
= V
.dimension()
976 for i
in xrange(1, V
.dimension()):
977 this_dim
= (u
**i
).operator_matrix().image().dimension()
978 if this_dim
< minimal_dim
:
979 minimal_dim
= this_dim
982 # Now minimal_matrix should correspond to the smallest
983 # non-zero subspace in Baes's (or really, Koecher's)
986 # However, we need to restrict the matrix to work on the
987 # subspace... or do we? Can't we just solve, knowing that
988 # A(c) = u^(s+1) should have a solution in the big space,
991 # Beware, solve_right() means that we're using COLUMN vectors.
992 # Our FiniteDimensionalAlgebraElement superclass uses rows.
994 A
= u_next
.operator_matrix()
995 c_coordinates
= A
.solve_right(u_next
.vector())
997 # Now c_coordinates is the idempotent we want, but it's in
998 # the coordinate system of the subalgebra.
1000 # We need the basis for J, but as elements of the parent algebra.
1002 basis
= [self
.parent(v
) for v
in V
.basis()]
1003 return self
.parent().linear_combination(zip(c_coordinates
, basis
))
1008 Return my trace, the sum of my eigenvalues.
1012 sage: J = JordanSpinEJA(3)
1013 sage: e0,e1,e2 = J.gens()
1014 sage: x = e0 + e1 + e2
1019 cs
= self
.characteristic_polynomial().coefficients(sparse
=False)
1023 raise ValueError('charpoly had fewer than 2 coefficients')
1026 def trace_inner_product(self
, other
):
1028 Return the trace inner product of myself and ``other``.
1030 if not other
in self
.parent():
1031 raise TypeError("'other' must live in the same algebra")
1033 return (self
*other
).trace()
1036 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1038 Return the Euclidean Jordan Algebra corresponding to the set
1039 `R^n` under the Hadamard product.
1041 Note: this is nothing more than the Cartesian product of ``n``
1042 copies of the spin algebra. Once Cartesian product algebras
1043 are implemented, this can go.
1047 This multiplication table can be verified by hand::
1049 sage: J = RealCartesianProductEJA(3)
1050 sage: e0,e1,e2 = J.gens()
1066 def __classcall_private__(cls
, n
, field
=QQ
):
1067 # The FiniteDimensionalAlgebra constructor takes a list of
1068 # matrices, the ith representing right multiplication by the ith
1069 # basis element in the vector space. So if e_1 = (1,0,0), then
1070 # right (Hadamard) multiplication of x by e_1 picks out the first
1071 # component of x; and likewise for the ith basis element e_i.
1072 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
1073 for i
in xrange(n
) ]
1075 fdeja
= super(RealCartesianProductEJA
, cls
)
1076 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=n
)
1078 def inner_product(self
, x
, y
):
1079 return _usual_ip(x
,y
)
1084 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1088 For now, we choose a random natural number ``n`` (greater than zero)
1089 and then give you back one of the following:
1091 * The cartesian product of the rational numbers ``n`` times; this is
1092 ``QQ^n`` with the Hadamard product.
1094 * The Jordan spin algebra on ``QQ^n``.
1096 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1099 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1100 in the space of ``2n``-by-``2n`` real symmetric matrices.
1102 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1103 in the space of ``4n``-by-``4n`` real symmetric matrices.
1105 Later this might be extended to return Cartesian products of the
1111 Euclidean Jordan algebra of degree...
1114 n
= ZZ
.random_element(1,5)
1115 constructor
= choice([RealCartesianProductEJA
,
1118 ComplexHermitianEJA
,
1119 QuaternionHermitianEJA
])
1120 return constructor(n
, field
=QQ
)
1124 def _real_symmetric_basis(n
, field
=QQ
):
1126 Return a basis for the space of real symmetric n-by-n matrices.
1128 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1132 for j
in xrange(i
+1):
1133 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1137 # Beware, orthogonal but not normalized!
1138 Sij
= Eij
+ Eij
.transpose()
1143 def _complex_hermitian_basis(n
, field
=QQ
):
1145 Returns a basis for the space of complex Hermitian n-by-n matrices.
1149 sage: set_random_seed()
1150 sage: n = ZZ.random_element(1,5)
1151 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1155 F
= QuadraticField(-1, 'I')
1158 # This is like the symmetric case, but we need to be careful:
1160 # * We want conjugate-symmetry, not just symmetry.
1161 # * The diagonal will (as a result) be real.
1165 for j
in xrange(i
+1):
1166 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1168 Sij
= _embed_complex_matrix(Eij
)
1171 # Beware, orthogonal but not normalized! The second one
1172 # has a minus because it's conjugated.
1173 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
1175 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
1180 def _quaternion_hermitian_basis(n
, field
=QQ
):
1182 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1186 sage: set_random_seed()
1187 sage: n = ZZ.random_element(1,5)
1188 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1192 Q
= QuaternionAlgebra(QQ
,-1,-1)
1195 # This is like the symmetric case, but we need to be careful:
1197 # * We want conjugate-symmetry, not just symmetry.
1198 # * The diagonal will (as a result) be real.
1202 for j
in xrange(i
+1):
1203 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1205 Sij
= _embed_quaternion_matrix(Eij
)
1208 # Beware, orthogonal but not normalized! The second,
1209 # third, and fourth ones have a minus because they're
1211 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
1213 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
1215 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
1217 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
1223 return vector(m
.base_ring(), m
.list())
1226 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
1228 def _multiplication_table_from_matrix_basis(basis
):
1230 At least three of the five simple Euclidean Jordan algebras have the
1231 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1232 multiplication on the right is matrix multiplication. Given a basis
1233 for the underlying matrix space, this function returns a
1234 multiplication table (obtained by looping through the basis
1235 elements) for an algebra of those matrices. A reordered copy
1236 of the basis is also returned to work around the fact that
1237 the ``span()`` in this function will change the order of the basis
1238 from what we think it is, to... something else.
1240 # In S^2, for example, we nominally have four coordinates even
1241 # though the space is of dimension three only. The vector space V
1242 # is supposed to hold the entire long vector, and the subspace W
1243 # of V will be spanned by the vectors that arise from symmetric
1244 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1245 field
= basis
[0].base_ring()
1246 dimension
= basis
[0].nrows()
1248 V
= VectorSpace(field
, dimension
**2)
1249 W
= V
.span( _mat2vec(s
) for s
in basis
)
1251 # Taking the span above reorders our basis (thanks, jerk!) so we
1252 # need to put our "matrix basis" in the same order as the
1253 # (reordered) vector basis.
1254 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
1258 # Brute force the multiplication-by-s matrix by looping
1259 # through all elements of the basis and doing the computation
1260 # to find out what the corresponding row should be. BEWARE:
1261 # these multiplication tables won't be symmetric! It therefore
1262 # becomes REALLY IMPORTANT that the underlying algebra
1263 # constructor uses ROW vectors and not COLUMN vectors. That's
1264 # why we're computing rows here and not columns.
1267 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
1268 Q_rows
.append(W
.coordinates(this_row
))
1269 Q
= matrix(field
, W
.dimension(), Q_rows
)
1275 def _embed_complex_matrix(M
):
1277 Embed the n-by-n complex matrix ``M`` into the space of real
1278 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1279 bi` to the block matrix ``[[a,b],[-b,a]]``.
1283 sage: F = QuadraticField(-1,'i')
1284 sage: x1 = F(4 - 2*i)
1285 sage: x2 = F(1 + 2*i)
1288 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1289 sage: _embed_complex_matrix(M)
1298 Embedding is a homomorphism (isomorphism, in fact)::
1300 sage: set_random_seed()
1301 sage: n = ZZ.random_element(5)
1302 sage: F = QuadraticField(-1, 'i')
1303 sage: X = random_matrix(F, n)
1304 sage: Y = random_matrix(F, n)
1305 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1306 sage: expected = _embed_complex_matrix(X*Y)
1307 sage: actual == expected
1313 raise ValueError("the matrix 'M' must be square")
1314 field
= M
.base_ring()
1319 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1321 # We can drop the imaginaries here.
1322 return block_matrix(field
.base_ring(), n
, blocks
)
1325 def _unembed_complex_matrix(M
):
1327 The inverse of _embed_complex_matrix().
1331 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1332 ....: [-2, 1, -4, 3],
1333 ....: [ 9, 10, 11, 12],
1334 ....: [-10, 9, -12, 11] ])
1335 sage: _unembed_complex_matrix(A)
1337 [ 10*i + 9 12*i + 11]
1341 Unembedding is the inverse of embedding::
1343 sage: set_random_seed()
1344 sage: F = QuadraticField(-1, 'i')
1345 sage: M = random_matrix(F, 3)
1346 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1352 raise ValueError("the matrix 'M' must be square")
1353 if not n
.mod(2).is_zero():
1354 raise ValueError("the matrix 'M' must be a complex embedding")
1356 F
= QuadraticField(-1, 'i')
1359 # Go top-left to bottom-right (reading order), converting every
1360 # 2-by-2 block we see to a single complex element.
1362 for k
in xrange(n
/2):
1363 for j
in xrange(n
/2):
1364 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1365 if submat
[0,0] != submat
[1,1]:
1366 raise ValueError('bad on-diagonal submatrix')
1367 if submat
[0,1] != -submat
[1,0]:
1368 raise ValueError('bad off-diagonal submatrix')
1369 z
= submat
[0,0] + submat
[0,1]*i
1372 return matrix(F
, n
/2, elements
)
1375 def _embed_quaternion_matrix(M
):
1377 Embed the n-by-n quaternion matrix ``M`` into the space of real
1378 matrices of size 4n-by-4n by first sending each quaternion entry
1379 `z = a + bi + cj + dk` to the block-complex matrix
1380 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1385 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1386 sage: i,j,k = Q.gens()
1387 sage: x = 1 + 2*i + 3*j + 4*k
1388 sage: M = matrix(Q, 1, [[x]])
1389 sage: _embed_quaternion_matrix(M)
1395 Embedding is a homomorphism (isomorphism, in fact)::
1397 sage: set_random_seed()
1398 sage: n = ZZ.random_element(5)
1399 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1400 sage: X = random_matrix(Q, n)
1401 sage: Y = random_matrix(Q, n)
1402 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1403 sage: expected = _embed_quaternion_matrix(X*Y)
1404 sage: actual == expected
1408 quaternions
= M
.base_ring()
1411 raise ValueError("the matrix 'M' must be square")
1413 F
= QuadraticField(-1, 'i')
1418 t
= z
.coefficient_tuple()
1423 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1424 [-c
+ d
*i
, a
- b
*i
]])
1425 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1427 # We should have real entries by now, so use the realest field
1428 # we've got for the return value.
1429 return block_matrix(quaternions
.base_ring(), n
, blocks
)
1432 def _unembed_quaternion_matrix(M
):
1434 The inverse of _embed_quaternion_matrix().
1438 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1439 ....: [-2, 1, -4, 3],
1440 ....: [-3, 4, 1, -2],
1441 ....: [-4, -3, 2, 1]])
1442 sage: _unembed_quaternion_matrix(M)
1443 [1 + 2*i + 3*j + 4*k]
1447 Unembedding is the inverse of embedding::
1449 sage: set_random_seed()
1450 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1451 sage: M = random_matrix(Q, 3)
1452 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1458 raise ValueError("the matrix 'M' must be square")
1459 if not n
.mod(4).is_zero():
1460 raise ValueError("the matrix 'M' must be a complex embedding")
1462 Q
= QuaternionAlgebra(QQ
,-1,-1)
1465 # Go top-left to bottom-right (reading order), converting every
1466 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1469 for l
in xrange(n
/4):
1470 for m
in xrange(n
/4):
1471 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1472 if submat
[0,0] != submat
[1,1].conjugate():
1473 raise ValueError('bad on-diagonal submatrix')
1474 if submat
[0,1] != -submat
[1,0].conjugate():
1475 raise ValueError('bad off-diagonal submatrix')
1476 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1477 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1480 return matrix(Q
, n
/4, elements
)
1483 # The usual inner product on R^n.
1485 return x
.vector().inner_product(y
.vector())
1487 # The inner product used for the real symmetric simple EJA.
1488 # We keep it as a separate function because e.g. the complex
1489 # algebra uses the same inner product, except divided by 2.
1490 def _matrix_ip(X
,Y
):
1491 X_mat
= X
.natural_representation()
1492 Y_mat
= Y
.natural_representation()
1493 return (X_mat
*Y_mat
).trace()
1496 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1498 The rank-n simple EJA consisting of real symmetric n-by-n
1499 matrices, the usual symmetric Jordan product, and the trace inner
1500 product. It has dimension `(n^2 + n)/2` over the reals.
1504 sage: J = RealSymmetricEJA(2)
1505 sage: e0, e1, e2 = J.gens()
1515 The degree of this algebra is `(n^2 + n) / 2`::
1517 sage: set_random_seed()
1518 sage: n = ZZ.random_element(1,5)
1519 sage: J = RealSymmetricEJA(n)
1520 sage: J.degree() == (n^2 + n)/2
1523 The Jordan multiplication is what we think it is::
1525 sage: set_random_seed()
1526 sage: n = ZZ.random_element(1,5)
1527 sage: J = RealSymmetricEJA(n)
1528 sage: x = J.random_element()
1529 sage: y = J.random_element()
1530 sage: actual = (x*y).natural_representation()
1531 sage: X = x.natural_representation()
1532 sage: Y = y.natural_representation()
1533 sage: expected = (X*Y + Y*X)/2
1534 sage: actual == expected
1536 sage: J(expected) == x*y
1541 def __classcall_private__(cls
, n
, field
=QQ
):
1542 S
= _real_symmetric_basis(n
, field
=field
)
1543 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
1545 fdeja
= super(RealSymmetricEJA
, cls
)
1546 return fdeja
.__classcall
_private
__(cls
,
1552 def inner_product(self
, x
, y
):
1553 return _matrix_ip(x
,y
)
1556 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1558 The rank-n simple EJA consisting of complex Hermitian n-by-n
1559 matrices over the real numbers, the usual symmetric Jordan product,
1560 and the real-part-of-trace inner product. It has dimension `n^2` over
1565 The degree of this algebra is `n^2`::
1567 sage: set_random_seed()
1568 sage: n = ZZ.random_element(1,5)
1569 sage: J = ComplexHermitianEJA(n)
1570 sage: J.degree() == n^2
1573 The Jordan multiplication is what we think it is::
1575 sage: set_random_seed()
1576 sage: n = ZZ.random_element(1,5)
1577 sage: J = ComplexHermitianEJA(n)
1578 sage: x = J.random_element()
1579 sage: y = J.random_element()
1580 sage: actual = (x*y).natural_representation()
1581 sage: X = x.natural_representation()
1582 sage: Y = y.natural_representation()
1583 sage: expected = (X*Y + Y*X)/2
1584 sage: actual == expected
1586 sage: J(expected) == x*y
1591 def __classcall_private__(cls
, n
, field
=QQ
):
1592 S
= _complex_hermitian_basis(n
)
1593 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
1595 fdeja
= super(ComplexHermitianEJA
, cls
)
1596 return fdeja
.__classcall
_private
__(cls
,
1602 def inner_product(self
, x
, y
):
1603 # Since a+bi on the diagonal is represented as
1608 # we'll double-count the "a" entries if we take the trace of
1610 return _matrix_ip(x
,y
)/2
1613 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1615 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1616 matrices, the usual symmetric Jordan product, and the
1617 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1622 The degree of this algebra is `n^2`::
1624 sage: set_random_seed()
1625 sage: n = ZZ.random_element(1,5)
1626 sage: J = QuaternionHermitianEJA(n)
1627 sage: J.degree() == 2*(n^2) - n
1630 The Jordan multiplication is what we think it is::
1632 sage: set_random_seed()
1633 sage: n = ZZ.random_element(1,5)
1634 sage: J = QuaternionHermitianEJA(n)
1635 sage: x = J.random_element()
1636 sage: y = J.random_element()
1637 sage: actual = (x*y).natural_representation()
1638 sage: X = x.natural_representation()
1639 sage: Y = y.natural_representation()
1640 sage: expected = (X*Y + Y*X)/2
1641 sage: actual == expected
1643 sage: J(expected) == x*y
1648 def __classcall_private__(cls
, n
, field
=QQ
):
1649 S
= _quaternion_hermitian_basis(n
)
1650 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
1652 fdeja
= super(QuaternionHermitianEJA
, cls
)
1653 return fdeja
.__classcall
_private
__(cls
,
1659 def inner_product(self
, x
, y
):
1660 # Since a+bi+cj+dk on the diagonal is represented as
1662 # a + bi +cj + dk = [ a b c d]
1667 # we'll quadruple-count the "a" entries if we take the trace of
1669 return _matrix_ip(x
,y
)/4
1672 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1674 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1675 with the usual inner product and jordan product ``x*y =
1676 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1681 This multiplication table can be verified by hand::
1683 sage: J = JordanSpinEJA(4)
1684 sage: e0,e1,e2,e3 = J.gens()
1702 def __classcall_private__(cls
, n
, field
=QQ
):
1704 id_matrix
= identity_matrix(field
, n
)
1706 ei
= id_matrix
.column(i
)
1707 Qi
= zero_matrix(field
, n
)
1709 Qi
.set_column(0, ei
)
1710 Qi
+= diagonal_matrix(n
, [ei
[0]]*n
)
1711 # The addition of the diagonal matrix adds an extra ei[0] in the
1712 # upper-left corner of the matrix.
1713 Qi
[0,0] = Qi
[0,0] * ~
field(2)
1716 # The rank of the spin algebra is two, unless we're in a
1717 # one-dimensional ambient space (because the rank is bounded by
1718 # the ambient dimension).
1719 fdeja
= super(JordanSpinEJA
, cls
)
1720 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=min(n
,2))
1722 def inner_product(self
, x
, y
):
1723 return _usual_ip(x
,y
)