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,
27 mult_table
= [b
.base_extend(field
) for b
in mult_table
]
30 if not (is_Matrix(b
) and b
.dimensions() == (n
, n
)):
31 raise ValueError("input is not a multiplication table")
32 mult_table
= tuple(mult_table
)
34 cat
= MagmaticAlgebras(field
).FiniteDimensional().WithBasis()
35 cat
.or_subcategory(category
)
36 if assume_associative
:
37 cat
= cat
.Associative()
39 names
= normalize_names(n
, names
)
41 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, cls
)
42 return fda
.__classcall
__(cls
,
45 assume_associative
=assume_associative
,
49 natural_basis
=natural_basis
,
50 inner_product
=inner_product
)
53 def __init__(self
, field
,
56 assume_associative
=False,
64 By definition, Jordan multiplication commutes::
66 sage: set_random_seed()
67 sage: J = random_eja()
68 sage: x = J.random_element()
69 sage: y = J.random_element()
75 self
._natural
_basis
= natural_basis
76 self
._inner
_product
= inner_product
77 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
86 Return a string representation of ``self``.
88 fmt
= "Euclidean Jordan algebra of degree {} over {}"
89 return fmt
.format(self
.degree(), self
.base_ring())
92 def inner_product(self
, x
, y
):
94 The inner product associated with this Euclidean Jordan algebra.
96 Will default to the trace inner product if nothing else.
98 if (not x
in self
) or (not y
in self
):
99 raise TypeError("arguments must live in this algebra")
100 if self
._inner
_product
is None:
101 return x
.trace_inner_product(y
)
103 return self
._inner
_product
(x
,y
)
106 def natural_basis(self
):
108 Return a more-natural representation of this algebra's basis.
110 Every finite-dimensional Euclidean Jordan Algebra is a direct
111 sum of five simple algebras, four of which comprise Hermitian
112 matrices. This method returns the original "natural" basis
113 for our underlying vector space. (Typically, the natural basis
114 is used to construct the multiplication table in the first place.)
116 Note that this will always return a matrix. The standard basis
117 in `R^n` will be returned as `n`-by-`1` column matrices.
121 sage: J = RealSymmetricSimpleEJA(2)
124 sage: J.natural_basis()
132 sage: J = JordanSpinSimpleEJA(2)
135 sage: J.natural_basis()
142 if self
._natural
_basis
is None:
143 return tuple( b
.vector().column() for b
in self
.basis() )
145 return self
._natural
_basis
150 Return the rank of this EJA.
152 if self
._rank
is None:
153 raise ValueError("no rank specified at genesis")
158 class Element(FiniteDimensionalAlgebraElement
):
160 An element of a Euclidean Jordan algebra.
163 def __init__(self
, A
, elt
=None):
167 The identity in `S^n` is converted to the identity in the EJA::
169 sage: J = RealSymmetricSimpleEJA(3)
170 sage: I = identity_matrix(QQ,3)
171 sage: J(I) == J.one()
174 This skew-symmetric matrix can't be represented in the EJA::
176 sage: J = RealSymmetricSimpleEJA(3)
177 sage: A = matrix(QQ,3, lambda i,j: i-j)
179 Traceback (most recent call last):
181 ArithmeticError: vector is not in free module
184 # Goal: if we're given a matrix, and if it lives in our
185 # parent algebra's "natural ambient space," convert it
186 # into an algebra element.
188 # The catch is, we make a recursive call after converting
189 # the given matrix into a vector that lives in the algebra.
190 # This we need to try the parent class initializer first,
191 # to avoid recursing forever if we're given something that
192 # already fits into the algebra, but also happens to live
193 # in the parent's "natural ambient space" (this happens with
196 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, elt
)
198 natural_basis
= A
.natural_basis()
199 if elt
in natural_basis
[0].matrix_space():
200 # Thanks for nothing! Matrix spaces aren't vector
201 # spaces in Sage, so we have to figure out its
202 # natural-basis coordinates ourselves.
203 V
= VectorSpace(elt
.base_ring(), elt
.nrows()**2)
204 W
= V
.span( _mat2vec(s
) for s
in natural_basis
)
205 coords
= W
.coordinates(_mat2vec(elt
))
206 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, coords
)
208 def __pow__(self
, n
):
210 Return ``self`` raised to the power ``n``.
212 Jordan algebras are always power-associative; see for
213 example Faraut and Koranyi, Proposition II.1.2 (ii).
217 We have to override this because our superclass uses row vectors
218 instead of column vectors! We, on the other hand, assume column
223 sage: set_random_seed()
224 sage: x = random_eja().random_element()
225 sage: x.operator_matrix()*x.vector() == (x^2).vector()
228 A few examples of power-associativity::
230 sage: set_random_seed()
231 sage: x = random_eja().random_element()
232 sage: x*(x*x)*(x*x) == x^5
234 sage: (x*x)*(x*x*x) == x^5
237 We also know that powers operator-commute (Koecher, Chapter
240 sage: set_random_seed()
241 sage: x = random_eja().random_element()
242 sage: m = ZZ.random_element(0,10)
243 sage: n = ZZ.random_element(0,10)
244 sage: Lxm = (x^m).operator_matrix()
245 sage: Lxn = (x^n).operator_matrix()
246 sage: Lxm*Lxn == Lxn*Lxm
256 return A( (self
.operator_matrix()**(n
-1))*self
.vector() )
259 def characteristic_polynomial(self
):
261 Return my characteristic polynomial (if I'm a regular
264 Eventually this should be implemented in terms of the parent
265 algebra's characteristic polynomial that works for ALL
268 if self
.is_regular():
269 return self
.minimal_polynomial()
271 raise NotImplementedError('irregular element')
274 def inner_product(self
, other
):
276 Return the parent algebra's inner product of myself and ``other``.
280 The inner product in the Jordan spin algebra is the usual
281 inner product on `R^n` (this example only works because the
282 basis for the Jordan algebra is the standard basis in `R^n`)::
284 sage: J = JordanSpinSimpleEJA(3)
285 sage: x = vector(QQ,[1,2,3])
286 sage: y = vector(QQ,[4,5,6])
287 sage: x.inner_product(y)
289 sage: J(x).inner_product(J(y))
292 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
293 multiplication is the usual matrix multiplication in `S^n`,
294 so the inner product of the identity matrix with itself
297 sage: J = RealSymmetricSimpleEJA(3)
298 sage: J.one().inner_product(J.one())
301 Likewise, the inner product on `C^n` is `<X,Y> =
302 Re(trace(X*Y))`, where we must necessarily take the real
303 part because the product of Hermitian matrices may not be
306 sage: J = ComplexHermitianSimpleEJA(3)
307 sage: J.one().inner_product(J.one())
312 Ensure that we can always compute an inner product, and that
313 it gives us back a real number::
315 sage: set_random_seed()
316 sage: J = random_eja()
317 sage: x = J.random_element()
318 sage: y = J.random_element()
319 sage: x.inner_product(y) in RR
325 raise TypeError("'other' must live in the same algebra")
327 return P
.inner_product(self
, other
)
330 def operator_commutes_with(self
, other
):
332 Return whether or not this element operator-commutes
337 The definition of a Jordan algebra says that any element
338 operator-commutes with its square::
340 sage: set_random_seed()
341 sage: x = random_eja().random_element()
342 sage: x.operator_commutes_with(x^2)
347 Test Lemma 1 from Chapter III of Koecher::
349 sage: set_random_seed()
350 sage: J = random_eja()
351 sage: u = J.random_element()
352 sage: v = J.random_element()
353 sage: lhs = u.operator_commutes_with(u*v)
354 sage: rhs = v.operator_commutes_with(u^2)
359 if not other
in self
.parent():
360 raise TypeError("'other' must live in the same algebra")
362 A
= self
.operator_matrix()
363 B
= other
.operator_matrix()
369 Return my determinant, the product of my eigenvalues.
373 sage: J = JordanSpinSimpleEJA(2)
374 sage: e0,e1 = J.gens()
378 sage: J = JordanSpinSimpleEJA(3)
379 sage: e0,e1,e2 = J.gens()
380 sage: x = e0 + e1 + e2
385 cs
= self
.characteristic_polynomial().coefficients(sparse
=False)
388 return cs
[0] * (-1)**r
390 raise ValueError('charpoly had no coefficients')
395 Return the Jordan-multiplicative inverse of this element.
397 We can't use the superclass method because it relies on the
398 algebra being associative.
402 The inverse in the spin factor algebra is given in Alizadeh's
405 sage: set_random_seed()
406 sage: n = ZZ.random_element(1,10)
407 sage: J = JordanSpinSimpleEJA(n)
408 sage: x = J.random_element()
409 sage: while x.is_zero():
410 ....: x = J.random_element()
411 sage: x_vec = x.vector()
413 sage: x_bar = x_vec[1:]
414 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
415 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
416 sage: x_inverse = coeff*inv_vec
417 sage: x.inverse() == J(x_inverse)
422 The identity element is its own inverse::
424 sage: set_random_seed()
425 sage: J = random_eja()
426 sage: J.one().inverse() == J.one()
429 If an element has an inverse, it acts like one. TODO: this
430 can be a lot less ugly once ``is_invertible`` doesn't crash
431 on irregular elements::
433 sage: set_random_seed()
434 sage: J = random_eja()
435 sage: x = J.random_element()
437 ....: x.inverse()*x == J.one()
443 if self
.parent().is_associative():
444 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
447 # TODO: we can do better once the call to is_invertible()
448 # doesn't crash on irregular elements.
449 #if not self.is_invertible():
450 # raise ValueError('element is not invertible')
452 # We do this a little different than the usual recursive
453 # call to a finite-dimensional algebra element, because we
454 # wind up with an inverse that lives in the subalgebra and
455 # we need information about the parent to convert it back.
456 V
= self
.span_of_powers()
457 assoc_subalg
= self
.subalgebra_generated_by()
458 # Mis-design warning: the basis used for span_of_powers()
459 # and subalgebra_generated_by() must be the same, and in
461 elt
= assoc_subalg(V
.coordinates(self
.vector()))
463 # This will be in the subalgebra's coordinates...
464 fda_elt
= FiniteDimensionalAlgebraElement(assoc_subalg
, elt
)
465 subalg_inverse
= fda_elt
.inverse()
467 # So we have to convert back...
468 basis
= [ self
.parent(v
) for v
in V
.basis() ]
469 pairs
= zip(subalg_inverse
.vector(), basis
)
470 return self
.parent().linear_combination(pairs
)
473 def is_invertible(self
):
475 Return whether or not this element is invertible.
477 We can't use the superclass method because it relies on
478 the algebra being associative.
480 return not self
.det().is_zero()
483 def is_nilpotent(self
):
485 Return whether or not some power of this element is zero.
487 The superclass method won't work unless we're in an
488 associative algebra, and we aren't. However, we generate
489 an assocoative subalgebra and we're nilpotent there if and
490 only if we're nilpotent here (probably).
494 The identity element is never nilpotent::
496 sage: set_random_seed()
497 sage: random_eja().one().is_nilpotent()
500 The additive identity is always nilpotent::
502 sage: set_random_seed()
503 sage: random_eja().zero().is_nilpotent()
507 # The element we're going to call "is_nilpotent()" on.
508 # Either myself, interpreted as an element of a finite-
509 # dimensional algebra, or an element of an associative
513 if self
.parent().is_associative():
514 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
516 V
= self
.span_of_powers()
517 assoc_subalg
= self
.subalgebra_generated_by()
518 # Mis-design warning: the basis used for span_of_powers()
519 # and subalgebra_generated_by() must be the same, and in
521 elt
= assoc_subalg(V
.coordinates(self
.vector()))
523 # Recursive call, but should work since elt lives in an
524 # associative algebra.
525 return elt
.is_nilpotent()
528 def is_regular(self
):
530 Return whether or not this is a regular element.
534 The identity element always has degree one, but any element
535 linearly-independent from it is regular::
537 sage: J = JordanSpinSimpleEJA(5)
538 sage: J.one().is_regular()
540 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
541 sage: for x in J.gens():
542 ....: (J.one() + x).is_regular()
550 return self
.degree() == self
.parent().rank()
555 Compute the degree of this element the straightforward way
556 according to the definition; by appending powers to a list
557 and figuring out its dimension (that is, whether or not
558 they're linearly dependent).
562 sage: J = JordanSpinSimpleEJA(4)
563 sage: J.one().degree()
565 sage: e0,e1,e2,e3 = J.gens()
566 sage: (e0 - e1).degree()
569 In the spin factor algebra (of rank two), all elements that
570 aren't multiples of the identity are regular::
572 sage: set_random_seed()
573 sage: n = ZZ.random_element(1,10)
574 sage: J = JordanSpinSimpleEJA(n)
575 sage: x = J.random_element()
576 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
580 return self
.span_of_powers().dimension()
583 def minimal_polynomial(self
):
587 sage: set_random_seed()
588 sage: x = random_eja().random_element()
589 sage: x.degree() == x.minimal_polynomial().degree()
594 sage: set_random_seed()
595 sage: x = random_eja().random_element()
596 sage: x.degree() == x.minimal_polynomial().degree()
599 The minimal polynomial and the characteristic polynomial coincide
600 and are known (see Alizadeh, Example 11.11) for all elements of
601 the spin factor algebra that aren't scalar multiples of the
604 sage: set_random_seed()
605 sage: n = ZZ.random_element(2,10)
606 sage: J = JordanSpinSimpleEJA(n)
607 sage: y = J.random_element()
608 sage: while y == y.coefficient(0)*J.one():
609 ....: y = J.random_element()
610 sage: y0 = y.vector()[0]
611 sage: y_bar = y.vector()[1:]
612 sage: actual = y.minimal_polynomial()
613 sage: x = SR.symbol('x', domain='real')
614 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
615 sage: bool(actual == expected)
619 # The element we're going to call "minimal_polynomial()" on.
620 # Either myself, interpreted as an element of a finite-
621 # dimensional algebra, or an element of an associative
625 if self
.parent().is_associative():
626 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
628 V
= self
.span_of_powers()
629 assoc_subalg
= self
.subalgebra_generated_by()
630 # Mis-design warning: the basis used for span_of_powers()
631 # and subalgebra_generated_by() must be the same, and in
633 elt
= assoc_subalg(V
.coordinates(self
.vector()))
635 # Recursive call, but should work since elt lives in an
636 # associative algebra.
637 return elt
.minimal_polynomial()
640 def natural_representation(self
):
642 Return a more-natural representation of this element.
644 Every finite-dimensional Euclidean Jordan Algebra is a
645 direct sum of five simple algebras, four of which comprise
646 Hermitian matrices. This method returns the original
647 "natural" representation of this element as a Hermitian
648 matrix, if it has one. If not, you get the usual representation.
652 sage: J = ComplexHermitianSimpleEJA(3)
655 sage: J.one().natural_representation()
664 B
= self
.parent().natural_basis()
665 W
= B
[0].matrix_space()
666 return W
.linear_combination(zip(self
.vector(), B
))
669 def operator_matrix(self
):
671 Return the matrix that represents left- (or right-)
672 multiplication by this element in the parent algebra.
674 We have to override this because the superclass method
675 returns a matrix that acts on row vectors (that is, on
680 Test the first polarization identity from my notes, Koecher Chapter
681 III, or from Baes (2.3)::
683 sage: set_random_seed()
684 sage: J = random_eja()
685 sage: x = J.random_element()
686 sage: y = J.random_element()
687 sage: Lx = x.operator_matrix()
688 sage: Ly = y.operator_matrix()
689 sage: Lxx = (x*x).operator_matrix()
690 sage: Lxy = (x*y).operator_matrix()
691 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
694 Test the second polarization identity from my notes or from
697 sage: set_random_seed()
698 sage: J = random_eja()
699 sage: x = J.random_element()
700 sage: y = J.random_element()
701 sage: z = J.random_element()
702 sage: Lx = x.operator_matrix()
703 sage: Ly = y.operator_matrix()
704 sage: Lz = z.operator_matrix()
705 sage: Lzy = (z*y).operator_matrix()
706 sage: Lxy = (x*y).operator_matrix()
707 sage: Lxz = (x*z).operator_matrix()
708 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
711 Test the third polarization identity from my notes or from
714 sage: set_random_seed()
715 sage: J = random_eja()
716 sage: u = J.random_element()
717 sage: y = J.random_element()
718 sage: z = J.random_element()
719 sage: Lu = u.operator_matrix()
720 sage: Ly = y.operator_matrix()
721 sage: Lz = z.operator_matrix()
722 sage: Lzy = (z*y).operator_matrix()
723 sage: Luy = (u*y).operator_matrix()
724 sage: Luz = (u*z).operator_matrix()
725 sage: Luyz = (u*(y*z)).operator_matrix()
726 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
727 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
728 sage: bool(lhs == rhs)
732 fda_elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
733 return fda_elt
.matrix().transpose()
736 def quadratic_representation(self
, other
=None):
738 Return the quadratic representation of this element.
742 The explicit form in the spin factor algebra is given by
743 Alizadeh's Example 11.12::
745 sage: set_random_seed()
746 sage: n = ZZ.random_element(1,10)
747 sage: J = JordanSpinSimpleEJA(n)
748 sage: x = J.random_element()
749 sage: x_vec = x.vector()
751 sage: x_bar = x_vec[1:]
752 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
753 sage: B = 2*x0*x_bar.row()
754 sage: C = 2*x0*x_bar.column()
755 sage: D = identity_matrix(QQ, n-1)
756 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
757 sage: D = D + 2*x_bar.tensor_product(x_bar)
758 sage: Q = block_matrix(2,2,[A,B,C,D])
759 sage: Q == x.quadratic_representation()
762 Test all of the properties from Theorem 11.2 in Alizadeh::
764 sage: set_random_seed()
765 sage: J = random_eja()
766 sage: x = J.random_element()
767 sage: y = J.random_element()
771 sage: actual = x.quadratic_representation(y)
772 sage: expected = ( (x+y).quadratic_representation()
773 ....: -x.quadratic_representation()
774 ....: -y.quadratic_representation() ) / 2
775 sage: actual == expected
780 sage: alpha = QQ.random_element()
781 sage: actual = (alpha*x).quadratic_representation()
782 sage: expected = (alpha^2)*x.quadratic_representation()
783 sage: actual == expected
788 sage: Qy = y.quadratic_representation()
789 sage: actual = J(Qy*x.vector()).quadratic_representation()
790 sage: expected = Qy*x.quadratic_representation()*Qy
791 sage: actual == expected
796 sage: k = ZZ.random_element(1,10)
797 sage: actual = (x^k).quadratic_representation()
798 sage: expected = (x.quadratic_representation())^k
799 sage: actual == expected
805 elif not other
in self
.parent():
806 raise TypeError("'other' must live in the same algebra")
808 L
= self
.operator_matrix()
809 M
= other
.operator_matrix()
810 return ( L
*M
+ M
*L
- (self
*other
).operator_matrix() )
813 def span_of_powers(self
):
815 Return the vector space spanned by successive powers of
818 # The dimension of the subalgebra can't be greater than
819 # the big algebra, so just put everything into a list
820 # and let span() get rid of the excess.
821 V
= self
.vector().parent()
822 return V
.span( (self
**d
).vector() for d
in xrange(V
.dimension()) )
825 def subalgebra_generated_by(self
):
827 Return the associative subalgebra of the parent EJA generated
832 sage: set_random_seed()
833 sage: x = random_eja().random_element()
834 sage: x.subalgebra_generated_by().is_associative()
837 Squaring in the subalgebra should be the same thing as
838 squaring in the superalgebra::
840 sage: set_random_seed()
841 sage: x = random_eja().random_element()
842 sage: u = x.subalgebra_generated_by().random_element()
843 sage: u.operator_matrix()*u.vector() == (u**2).vector()
847 # First get the subspace spanned by the powers of myself...
848 V
= self
.span_of_powers()
851 # Now figure out the entries of the right-multiplication
852 # matrix for the successive basis elements b0, b1,... of
855 for b_right
in V
.basis():
856 eja_b_right
= self
.parent()(b_right
)
858 # The first row of the right-multiplication matrix by
859 # b1 is what we get if we apply that matrix to b1. The
860 # second row of the right multiplication matrix by b1
861 # is what we get when we apply that matrix to b2...
863 # IMPORTANT: this assumes that all vectors are COLUMN
864 # vectors, unlike our superclass (which uses row vectors).
865 for b_left
in V
.basis():
866 eja_b_left
= self
.parent()(b_left
)
867 # Multiply in the original EJA, but then get the
868 # coordinates from the subalgebra in terms of its
870 this_row
= V
.coordinates((eja_b_left
*eja_b_right
).vector())
871 b_right_rows
.append(this_row
)
872 b_right_matrix
= matrix(F
, b_right_rows
)
873 mats
.append(b_right_matrix
)
875 # It's an algebra of polynomials in one element, and EJAs
876 # are power-associative.
878 # TODO: choose generator names intelligently.
879 return FiniteDimensionalEuclideanJordanAlgebra(F
, mats
, assume_associative
=True, names
='f')
882 def subalgebra_idempotent(self
):
884 Find an idempotent in the associative subalgebra I generate
885 using Proposition 2.3.5 in Baes.
889 sage: set_random_seed()
891 sage: c = J.random_element().subalgebra_idempotent()
894 sage: J = JordanSpinSimpleEJA(5)
895 sage: c = J.random_element().subalgebra_idempotent()
900 if self
.is_nilpotent():
901 raise ValueError("this only works with non-nilpotent elements!")
903 V
= self
.span_of_powers()
904 J
= self
.subalgebra_generated_by()
905 # Mis-design warning: the basis used for span_of_powers()
906 # and subalgebra_generated_by() must be the same, and in
908 u
= J(V
.coordinates(self
.vector()))
910 # The image of the matrix of left-u^m-multiplication
911 # will be minimal for some natural number s...
913 minimal_dim
= V
.dimension()
914 for i
in xrange(1, V
.dimension()):
915 this_dim
= (u
**i
).operator_matrix().image().dimension()
916 if this_dim
< minimal_dim
:
917 minimal_dim
= this_dim
920 # Now minimal_matrix should correspond to the smallest
921 # non-zero subspace in Baes's (or really, Koecher's)
924 # However, we need to restrict the matrix to work on the
925 # subspace... or do we? Can't we just solve, knowing that
926 # A(c) = u^(s+1) should have a solution in the big space,
929 # Beware, solve_right() means that we're using COLUMN vectors.
930 # Our FiniteDimensionalAlgebraElement superclass uses rows.
932 A
= u_next
.operator_matrix()
933 c_coordinates
= A
.solve_right(u_next
.vector())
935 # Now c_coordinates is the idempotent we want, but it's in
936 # the coordinate system of the subalgebra.
938 # We need the basis for J, but as elements of the parent algebra.
940 basis
= [self
.parent(v
) for v
in V
.basis()]
941 return self
.parent().linear_combination(zip(c_coordinates
, basis
))
946 Return my trace, the sum of my eigenvalues.
950 sage: J = JordanSpinSimpleEJA(3)
951 sage: e0,e1,e2 = J.gens()
952 sage: x = e0 + e1 + e2
957 cs
= self
.characteristic_polynomial().coefficients(sparse
=False)
961 raise ValueError('charpoly had fewer than 2 coefficients')
964 def trace_inner_product(self
, other
):
966 Return the trace inner product of myself and ``other``.
968 if not other
in self
.parent():
969 raise TypeError("'other' must live in the same algebra")
971 return (self
*other
).trace()
974 def eja_rn(dimension
, field
=QQ
):
976 Return the Euclidean Jordan Algebra corresponding to the set
977 `R^n` under the Hadamard product.
981 This multiplication table can be verified by hand::
984 sage: e0,e1,e2 = J.gens()
999 # The FiniteDimensionalAlgebra constructor takes a list of
1000 # matrices, the ith representing right multiplication by the ith
1001 # basis element in the vector space. So if e_1 = (1,0,0), then
1002 # right (Hadamard) multiplication of x by e_1 picks out the first
1003 # component of x; and likewise for the ith basis element e_i.
1004 Qs
= [ matrix(field
, dimension
, dimension
, lambda k
,j
: 1*(k
== j
== i
))
1005 for i
in xrange(dimension
) ]
1007 return FiniteDimensionalEuclideanJordanAlgebra(field
,
1010 inner_product
=_usual_ip
)
1016 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1020 For now, we choose a random natural number ``n`` (greater than zero)
1021 and then give you back one of the following:
1023 * The cartesian product of the rational numbers ``n`` times; this is
1024 ``QQ^n`` with the Hadamard product.
1026 * The Jordan spin algebra on ``QQ^n``.
1028 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1031 Later this might be extended to return Cartesian products of the
1037 Euclidean Jordan algebra of degree...
1040 n
= ZZ
.random_element(1,5)
1041 constructor
= choice([eja_rn
,
1042 JordanSpinSimpleEJA
,
1043 RealSymmetricSimpleEJA
,
1044 ComplexHermitianSimpleEJA
])
1045 return constructor(n
, field
=QQ
)
1049 def _real_symmetric_basis(n
, field
=QQ
):
1051 Return a basis for the space of real symmetric n-by-n matrices.
1053 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1057 for j
in xrange(i
+1):
1058 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1062 # Beware, orthogonal but not normalized!
1063 Sij
= Eij
+ Eij
.transpose()
1068 def _complex_hermitian_basis(n
, field
=QQ
):
1070 Returns a basis for the space of complex Hermitian n-by-n matrices.
1074 sage: set_random_seed()
1075 sage: n = ZZ.random_element(1,5)
1076 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1080 F
= QuadraticField(-1, 'I')
1083 # This is like the symmetric case, but we need to be careful:
1085 # * We want conjugate-symmetry, not just symmetry.
1086 # * The diagonal will (as a result) be real.
1090 for j
in xrange(i
+1):
1091 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1093 Sij
= _embed_complex_matrix(Eij
)
1096 # Beware, orthogonal but not normalized! The second one
1097 # has a minus because it's conjugated.
1098 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
1100 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
1106 return vector(m
.base_ring(), m
.list())
1109 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
1111 def _multiplication_table_from_matrix_basis(basis
):
1113 At least three of the five simple Euclidean Jordan algebras have the
1114 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1115 multiplication on the right is matrix multiplication. Given a basis
1116 for the underlying matrix space, this function returns a
1117 multiplication table (obtained by looping through the basis
1118 elements) for an algebra of those matrices. A reordered copy
1119 of the basis is also returned to work around the fact that
1120 the ``span()`` in this function will change the order of the basis
1121 from what we think it is, to... something else.
1123 # In S^2, for example, we nominally have four coordinates even
1124 # though the space is of dimension three only. The vector space V
1125 # is supposed to hold the entire long vector, and the subspace W
1126 # of V will be spanned by the vectors that arise from symmetric
1127 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1128 field
= basis
[0].base_ring()
1129 dimension
= basis
[0].nrows()
1131 V
= VectorSpace(field
, dimension
**2)
1132 W
= V
.span( _mat2vec(s
) for s
in basis
)
1134 # Taking the span above reorders our basis (thanks, jerk!) so we
1135 # need to put our "matrix basis" in the same order as the
1136 # (reordered) vector basis.
1137 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
1141 # Brute force the multiplication-by-s matrix by looping
1142 # through all elements of the basis and doing the computation
1143 # to find out what the corresponding row should be. BEWARE:
1144 # these multiplication tables won't be symmetric! It therefore
1145 # becomes REALLY IMPORTANT that the underlying algebra
1146 # constructor uses ROW vectors and not COLUMN vectors. That's
1147 # why we're computing rows here and not columns.
1150 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
1151 Q_rows
.append(W
.coordinates(this_row
))
1152 Q
= matrix(field
, W
.dimension(), Q_rows
)
1158 def _embed_complex_matrix(M
):
1160 Embed the n-by-n complex matrix ``M`` into the space of real
1161 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1162 bi` to the block matrix ``[[a,b],[-b,a]]``.
1166 sage: F = QuadraticField(-1,'i')
1167 sage: x1 = F(4 - 2*i)
1168 sage: x2 = F(1 + 2*i)
1171 sage: M = matrix(F,2,[x1,x2,x3,x4])
1172 sage: _embed_complex_matrix(M)
1182 raise ValueError("the matrix 'M' must be square")
1183 field
= M
.base_ring()
1188 blocks
.append(matrix(field
, 2, [[a
,-b
],[b
,a
]]))
1190 # We can drop the imaginaries here.
1191 return block_matrix(field
.base_ring(), n
, blocks
)
1194 def _unembed_complex_matrix(M
):
1196 The inverse of _embed_complex_matrix().
1200 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1201 ....: [-2, 1, -4, 3],
1202 ....: [ 9, 10, 11, 12],
1203 ....: [-10, 9, -12, 11] ])
1204 sage: _unembed_complex_matrix(A)
1205 [ -2*i + 1 -4*i + 3]
1206 [ -10*i + 9 -12*i + 11]
1210 raise ValueError("the matrix 'M' must be square")
1211 if not n
.mod(2).is_zero():
1212 raise ValueError("the matrix 'M' must be a complex embedding")
1214 F
= QuadraticField(-1, 'i')
1217 # Go top-left to bottom-right (reading order), converting every
1218 # 2-by-2 block we see to a single complex element.
1220 for k
in xrange(n
/2):
1221 for j
in xrange(n
/2):
1222 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1223 if submat
[0,0] != submat
[1,1]:
1224 raise ValueError('bad real submatrix')
1225 if submat
[0,1] != -submat
[1,0]:
1226 raise ValueError('bad imag submatrix')
1227 z
= submat
[0,0] + submat
[1,0]*i
1230 return matrix(F
, n
/2, elements
)
1232 # The usual inner product on R^n.
1234 return x
.vector().inner_product(y
.vector())
1236 # The inner product used for the real symmetric simple EJA.
1237 # We keep it as a separate function because e.g. the complex
1238 # algebra uses the same inner product, except divided by 2.
1239 def _matrix_ip(X
,Y
):
1240 X_mat
= X
.natural_representation()
1241 Y_mat
= Y
.natural_representation()
1242 return (X_mat
*Y_mat
).trace()
1245 def RealSymmetricSimpleEJA(n
, field
=QQ
):
1247 The rank-n simple EJA consisting of real symmetric n-by-n
1248 matrices, the usual symmetric Jordan product, and the trace inner
1249 product. It has dimension `(n^2 + n)/2` over the reals.
1253 sage: J = RealSymmetricSimpleEJA(2)
1254 sage: e0, e1, e2 = J.gens()
1264 The degree of this algebra is `(n^2 + n) / 2`::
1266 sage: set_random_seed()
1267 sage: n = ZZ.random_element(1,5)
1268 sage: J = RealSymmetricSimpleEJA(n)
1269 sage: J.degree() == (n^2 + n)/2
1273 S
= _real_symmetric_basis(n
, field
=field
)
1274 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
1276 return FiniteDimensionalEuclideanJordanAlgebra(field
,
1280 inner_product
=_matrix_ip
)
1283 def ComplexHermitianSimpleEJA(n
, field
=QQ
):
1285 The rank-n simple EJA consisting of complex Hermitian n-by-n
1286 matrices over the real numbers, the usual symmetric Jordan product,
1287 and the real-part-of-trace inner product. It has dimension `n^2` over
1292 The degree of this algebra is `n^2`::
1294 sage: set_random_seed()
1295 sage: n = ZZ.random_element(1,5)
1296 sage: J = ComplexHermitianSimpleEJA(n)
1297 sage: J.degree() == n^2
1301 S
= _complex_hermitian_basis(n
)
1302 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
1304 # Since a+bi on the diagonal is represented as
1309 # we'll double-count the "a" entries if we take the trace of
1311 ip
= lambda X
,Y
: _matrix_ip(X
,Y
)/2
1313 return FiniteDimensionalEuclideanJordanAlgebra(field
,
1320 def QuaternionHermitianSimpleEJA(n
):
1322 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1323 matrices, the usual symmetric Jordan product, and the
1324 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1329 def OctonionHermitianSimpleEJA(n
):
1331 This shit be crazy. It has dimension 27 over the reals.
1336 def JordanSpinSimpleEJA(n
, field
=QQ
):
1338 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1339 with the usual inner product and jordan product ``x*y =
1340 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1345 This multiplication table can be verified by hand::
1347 sage: J = JordanSpinSimpleEJA(4)
1348 sage: e0,e1,e2,e3 = J.gens()
1364 In one dimension, this is the reals under multiplication::
1366 sage: J1 = JordanSpinSimpleEJA(1)
1367 sage: J2 = eja_rn(1)
1373 id_matrix
= identity_matrix(field
, n
)
1375 ei
= id_matrix
.column(i
)
1376 Qi
= zero_matrix(field
, n
)
1378 Qi
.set_column(0, ei
)
1379 Qi
+= diagonal_matrix(n
, [ei
[0]]*n
)
1380 # The addition of the diagonal matrix adds an extra ei[0] in the
1381 # upper-left corner of the matrix.
1382 Qi
[0,0] = Qi
[0,0] * ~
field(2)
1385 # The rank of the spin factor algebra is two, UNLESS we're in a
1386 # one-dimensional ambient space (the rank is bounded by the
1387 # ambient dimension).
1388 return FiniteDimensionalEuclideanJordanAlgebra(field
,
1391 inner_product
=_usual_ip
)