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
.finite_dimensional_algebras_with_basis
import FiniteDimensionalAlgebrasWithBasis
9 from sage
.categories
.morphism
import SetMorphism
10 from sage
.structure
.element
import is_Matrix
11 from sage
.structure
.category_object
import normalize_names
13 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra
import FiniteDimensionalAlgebra
14 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra_element
import FiniteDimensionalAlgebraElement
15 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra_morphism
import FiniteDimensionalAlgebraMorphism
, FiniteDimensionalAlgebraHomset
18 class FiniteDimensionalEuclideanJordanAlgebraHomset(FiniteDimensionalAlgebraHomset
):
20 def has_coerce_map_from(self
, S
):
24 sage: J = RealSymmetricEJA(2)
26 sage: H.has_coerce_map_from(QQ)
31 # The Homset classes override has_coerce_map_from() with
32 # something that crashes when it's given e.g. QQ.
33 if S
.is_subring(self
.codomain().base_ring()):
36 pclass
= super(FiniteDimensionalEuclideanJordanAlgebraHomset
, self
)
37 return pclass
.has_coerce_map_from(S
)
40 def _coerce_map_from_(self
, S
):
44 sage: J = RealSymmetricEJA(2)
47 Morphism from Euclidean Jordan algebra of degree 3 over Rational
48 Field to Euclidean Jordan algebra of degree 3 over Rational Field
58 h
= S
.hom(self
.codomain())
59 return SetMorphism(Hom(S
,C
), lambda x
: h(x
).operator())
62 def __call__(self
, x
):
66 sage: J = RealSymmetricEJA(2)
69 Morphism from Euclidean Jordan algebra of degree 3 over Rational
70 Field to Euclidean Jordan algebra of degree 3 over Rational Field
77 if x
in self
.base_ring():
78 cols
= self
.domain().dimension()
79 rows
= self
.codomain().dimension()
80 x
= x
*identity_matrix(self
.codomain().base_ring(), rows
, cols
)
81 return FiniteDimensionalEuclideanJordanAlgebraMorphism(self
, x
)
84 class FiniteDimensionalEuclideanJordanAlgebraMorphism(FiniteDimensionalAlgebraMorphism
):
86 A linear map between two finite-dimensional EJAs.
88 This is a very thin wrapper around FiniteDimensionalAlgebraMorphism
89 that does only a few things:
91 1. Avoids the ``unitary`` and ``check`` arguments to the constructor
92 that will always be ``False``. This is necessary because these
93 are homomorphisms with respect to ADDITION, but the SageMath
94 machinery wants to check that they're homomorphisms with respect
95 to (Jordan) MULTIPLICATION. That obviously doesn't work.
97 2. Inputs and outputs the underlying matrix with respect to COLUMN
98 vectors, unlike the parent class.
100 3. Allows us to add, subtract, negate, multiply (compose), and
101 invert morphisms in the obvious way.
103 If this seems a bit heavyweight, it is. I would have been happy to
104 use a the ring morphism that underlies the finite-dimensional
105 algebra morphism, but they don't seem to be callable on elements of
106 our EJA, and you can't add/multiply/etc. them.
108 def _add_(self
, other
):
110 Add two EJA morphisms in the obvious way.
114 sage: J = RealSymmetricEJA(3)
117 sage: x.operator() + y.operator()
118 Morphism from Euclidean Jordan algebra of degree 6 over Rational
119 Field to Euclidean Jordan algebra of degree 6 over Rational Field
130 sage: set_random_seed()
131 sage: J = random_eja()
132 sage: x = J.random_element()
133 sage: y = J.random_element()
134 sage: (x.operator() + y.operator()) in J.Hom(J)
140 raise ValueError("summands must live in the same space")
142 return FiniteDimensionalEuclideanJordanAlgebraMorphism(
144 self
.matrix() + other
.matrix() )
147 def __init__(self
, parent
, f
):
148 FiniteDimensionalAlgebraMorphism
.__init
__(self
,
155 def __invert__(self
):
159 sage: J = RealSymmetricEJA(2)
160 sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
161 sage: x.is_invertible()
164 Morphism from Euclidean Jordan algebra of degree 3 over Rational
165 Field to Euclidean Jordan algebra of degree 3 over Rational Field
170 sage: x.operator_matrix().inverse()
177 sage: set_random_seed()
178 sage: J = random_eja()
179 sage: x = J.random_element()
180 sage: not x.is_invertible() or (
181 ....: (~x.operator()).matrix() == x.operator_matrix().inverse() )
186 if not A
.is_invertible():
187 raise ValueError("morphism is not invertible")
190 return FiniteDimensionalEuclideanJordanAlgebraMorphism(self
.parent(),
193 def _lmul_(self
, right
):
195 Compose two EJA morphisms using multiplicative notation.
199 sage: J = RealSymmetricEJA(2)
202 sage: x.operator() * y.operator()
203 Morphism from Euclidean Jordan algebra of degree 3 over Rational
204 Field to Euclidean Jordan algebra of degree 3 over Rational Field
212 sage: J = RealSymmetricEJA(2)
213 sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
215 Morphism from Euclidean Jordan algebra of degree 3 over Rational
216 Field to Euclidean Jordan algebra of degree 3 over Rational Field
222 Morphism from Euclidean Jordan algebra of degree 3 over Rational
223 Field to Euclidean Jordan algebra of degree 3 over Rational Field
229 Morphism from Euclidean Jordan algebra of degree 3 over Rational
230 Field to Euclidean Jordan algebra of degree 3 over Rational Field
238 sage: set_random_seed()
239 sage: J = random_eja()
240 sage: x = J.random_element()
241 sage: y = J.random_element()
242 sage: (x.operator() * y.operator()) in J.Hom(J)
247 # I think the morphism classes break the coercion framework
248 # somewhere along the way, so we have to do this ourselves.
249 right
= self
.parent().coerce(right
)
253 if not right
.codomain() is self
.domain():
254 raise ValueError("(co)domains must agree for composition")
256 return FiniteDimensionalEuclideanJordanAlgebraMorphism(
258 self
.matrix()*right
.matrix() )
265 Negate this morphism.
269 sage: J = RealSymmetricEJA(2)
272 Morphism from Euclidean Jordan algebra of degree 3 over Rational
273 Field to Euclidean Jordan algebra of degree 3 over Rational Field
281 sage: set_random_seed()
282 sage: J = random_eja()
283 sage: x = J.random_element()
284 sage: -x.operator() in J.Hom(J)
288 return FiniteDimensionalEuclideanJordanAlgebraMorphism(
295 We override only the representation that is shown to the user,
296 because we want the matrix to be with respect to COLUMN vectors.
300 Ensure that we see the transpose of the underlying matrix object:
302 sage: J = RealSymmetricEJA(3)
303 sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
304 sage: L = x.operator()
306 Morphism from Euclidean Jordan algebra of degree 6 over Rational
307 Field to Euclidean Jordan algebra of degree 6 over Rational Field
324 return "Morphism from {} to {} given by matrix\n{}".format(
325 self
.domain(), self
.codomain(), self
.matrix())
328 def __sub__(self
, other
):
330 Subtract one morphism from another using addition and negation.
334 sage: J = RealSymmetricEJA(2)
335 sage: L1 = J.one().operator()
337 Morphism from Euclidean Jordan algebra of degree 3 over Rational
338 Field to Euclidean Jordan algebra of degree 3 over Rational
339 Field given by matrix
346 sage: set_random_seed()
347 sage: J = random_eja()
348 sage: x = J.random_element()
349 sage: y = J.random_element()
350 sage: x.operator() - y.operator() in J.Hom(J)
354 return self
+ (-other
)
359 Return the matrix of this morphism with respect to a left-action
362 return FiniteDimensionalAlgebraMorphism
.matrix(self
).transpose()
365 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra
):
367 def __classcall_private__(cls
,
371 assume_associative
=False,
376 mult_table
= [b
.base_extend(field
) for b
in mult_table
]
379 if not (is_Matrix(b
) and b
.dimensions() == (n
, n
)):
380 raise ValueError("input is not a multiplication table")
381 mult_table
= tuple(mult_table
)
383 cat
= FiniteDimensionalAlgebrasWithBasis(field
)
384 cat
.or_subcategory(category
)
385 if assume_associative
:
386 cat
= cat
.Associative()
388 names
= normalize_names(n
, names
)
390 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, cls
)
391 return fda
.__classcall
__(cls
,
394 assume_associative
=assume_associative
,
398 natural_basis
=natural_basis
)
401 def _Hom_(self
, B
, cat
):
403 Construct a homset of ``self`` and ``B``.
405 return FiniteDimensionalEuclideanJordanAlgebraHomset(self
,
414 assume_associative
=False,
421 By definition, Jordan multiplication commutes::
423 sage: set_random_seed()
424 sage: J = random_eja()
425 sage: x = J.random_element()
426 sage: y = J.random_element()
432 self
._natural
_basis
= natural_basis
433 self
._multiplication
_table
= mult_table
434 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
443 Return a string representation of ``self``.
445 fmt
= "Euclidean Jordan algebra of degree {} over {}"
446 return fmt
.format(self
.degree(), self
.base_ring())
449 def _a_regular_element(self
):
451 Guess a regular element. Needed to compute the basis for our
452 characteristic polynomial coefficients.
455 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
456 if not z
.is_regular():
457 raise ValueError("don't know a regular element")
462 def _charpoly_basis_space(self
):
464 Return the vector space spanned by the basis used in our
465 characteristic polynomial coefficients. This is used not only to
466 compute those coefficients, but also any time we need to
467 evaluate the coefficients (like when we compute the trace or
470 z
= self
._a
_regular
_element
()
471 V
= z
.vector().parent().ambient_vector_space()
472 V1
= V
.span_of_basis( (z
**k
).vector() for k
in range(self
.rank()) )
473 b
= (V1
.basis() + V1
.complement().basis())
474 return V
.span_of_basis(b
)
478 def _charpoly_coeff(self
, i
):
480 Return the coefficient polynomial "a_{i}" of this algebra's
481 general characteristic polynomial.
483 Having this be a separate cached method lets us compute and
484 store the trace/determinant (a_{r-1} and a_{0} respectively)
485 separate from the entire characteristic polynomial.
487 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
488 R
= A_of_x
.base_ring()
490 # Guaranteed by theory
493 # Danger: the in-place modification is done for performance
494 # reasons (reconstructing a matrix with huge polynomial
495 # entries is slow), but I don't know how cached_method works,
496 # so it's highly possible that we're modifying some global
497 # list variable by reference, here. In other words, you
498 # probably shouldn't call this method twice on the same
499 # algebra, at the same time, in two threads
500 Ai_orig
= A_of_x
.column(i
)
501 A_of_x
.set_column(i
,xr
)
502 numerator
= A_of_x
.det()
503 A_of_x
.set_column(i
,Ai_orig
)
505 # We're relying on the theory here to ensure that each a_i is
506 # indeed back in R, and the added negative signs are to make
507 # the whole charpoly expression sum to zero.
508 return R(-numerator
/detA
)
512 def _charpoly_matrix_system(self
):
514 Compute the matrix whose entries A_ij are polynomials in
515 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
516 corresponding to `x^r` and the determinent of the matrix A =
517 [A_ij]. In other words, all of the fixed (cachable) data needed
518 to compute the coefficients of the characteristic polynomial.
523 # Construct a new algebra over a multivariate polynomial ring...
524 names
= ['X' + str(i
) for i
in range(1,n
+1)]
525 R
= PolynomialRing(self
.base_ring(), names
)
526 J
= FiniteDimensionalEuclideanJordanAlgebra(R
,
527 self
._multiplication
_table
,
530 idmat
= identity_matrix(J
.base_ring(), n
)
532 W
= self
._charpoly
_basis
_space
()
533 W
= W
.change_ring(R
.fraction_field())
535 # Starting with the standard coordinates x = (X1,X2,...,Xn)
536 # and then converting the entries to W-coordinates allows us
537 # to pass in the standard coordinates to the charpoly and get
538 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
541 # W.coordinates(x^2) eval'd at (standard z-coords)
545 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
547 # We want the middle equivalent thing in our matrix, but use
548 # the first equivalent thing instead so that we can pass in
549 # standard coordinates.
550 x
= J(vector(R
, R
.gens()))
551 l1
= [column_matrix(W
.coordinates((x
**k
).vector())) for k
in range(r
)]
552 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
553 A_of_x
= block_matrix(R
, 1, n
, (l1
+ l2
))
554 xr
= W
.coordinates((x
**r
).vector())
555 return (A_of_x
, x
, xr
, A_of_x
.det())
559 def characteristic_polynomial(self
):
564 This implementation doesn't guarantee that the polynomial
565 denominator in the coefficients is not identically zero, so
566 theoretically it could crash. The way that this is handled
567 in e.g. Faraut and Koranyi is to use a basis that guarantees
568 the denominator is non-zero. But, doing so requires knowledge
569 of at least one regular element, and we don't even know how
570 to do that. The trade-off is that, if we use the standard basis,
571 the resulting polynomial will accept the "usual" coordinates. In
572 other words, we don't have to do a change of basis before e.g.
573 computing the trace or determinant.
577 The characteristic polynomial in the spin algebra is given in
578 Alizadeh, Example 11.11::
580 sage: J = JordanSpinEJA(3)
581 sage: p = J.characteristic_polynomial(); p
582 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
583 sage: xvec = J.one().vector()
591 # The list of coefficient polynomials a_1, a_2, ..., a_n.
592 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
594 # We go to a bit of trouble here to reorder the
595 # indeterminates, so that it's easier to evaluate the
596 # characteristic polynomial at x's coordinates and get back
597 # something in terms of t, which is what we want.
599 S
= PolynomialRing(self
.base_ring(),'t')
601 S
= PolynomialRing(S
, R
.variable_names())
604 # Note: all entries past the rth should be zero. The
605 # coefficient of the highest power (x^r) is 1, but it doesn't
606 # appear in the solution vector which contains coefficients
607 # for the other powers (to make them sum to x^r).
609 a
[r
] = 1 # corresponds to x^r
611 # When the rank is equal to the dimension, trying to
612 # assign a[r] goes out-of-bounds.
613 a
.append(1) # corresponds to x^r
615 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
618 def inner_product(self
, x
, y
):
620 The inner product associated with this Euclidean Jordan algebra.
622 Defaults to the trace inner product, but can be overridden by
623 subclasses if they are sure that the necessary properties are
628 The inner product must satisfy its axiom for this algebra to truly
629 be a Euclidean Jordan Algebra::
631 sage: set_random_seed()
632 sage: J = random_eja()
633 sage: x = J.random_element()
634 sage: y = J.random_element()
635 sage: z = J.random_element()
636 sage: (x*y).inner_product(z) == y.inner_product(x*z)
640 if (not x
in self
) or (not y
in self
):
641 raise TypeError("arguments must live in this algebra")
642 return x
.trace_inner_product(y
)
645 def natural_basis(self
):
647 Return a more-natural representation of this algebra's basis.
649 Every finite-dimensional Euclidean Jordan Algebra is a direct
650 sum of five simple algebras, four of which comprise Hermitian
651 matrices. This method returns the original "natural" basis
652 for our underlying vector space. (Typically, the natural basis
653 is used to construct the multiplication table in the first place.)
655 Note that this will always return a matrix. The standard basis
656 in `R^n` will be returned as `n`-by-`1` column matrices.
660 sage: J = RealSymmetricEJA(2)
663 sage: J.natural_basis()
671 sage: J = JordanSpinEJA(2)
674 sage: J.natural_basis()
681 if self
._natural
_basis
is None:
682 return tuple( b
.vector().column() for b
in self
.basis() )
684 return self
._natural
_basis
689 Return the rank of this EJA.
691 if self
._rank
is None:
692 raise ValueError("no rank specified at genesis")
697 class Element(FiniteDimensionalAlgebraElement
):
699 An element of a Euclidean Jordan algebra.
704 Oh man, I should not be doing this. This hides the "disabled"
705 methods ``left_matrix`` and ``matrix`` from introspection;
706 in particular it removes them from tab-completion.
708 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
709 dir(self
.__class
__) )
712 def __init__(self
, A
, elt
=None):
716 The identity in `S^n` is converted to the identity in the EJA::
718 sage: J = RealSymmetricEJA(3)
719 sage: I = identity_matrix(QQ,3)
720 sage: J(I) == J.one()
723 This skew-symmetric matrix can't be represented in the EJA::
725 sage: J = RealSymmetricEJA(3)
726 sage: A = matrix(QQ,3, lambda i,j: i-j)
728 Traceback (most recent call last):
730 ArithmeticError: vector is not in free module
733 # Goal: if we're given a matrix, and if it lives in our
734 # parent algebra's "natural ambient space," convert it
735 # into an algebra element.
737 # The catch is, we make a recursive call after converting
738 # the given matrix into a vector that lives in the algebra.
739 # This we need to try the parent class initializer first,
740 # to avoid recursing forever if we're given something that
741 # already fits into the algebra, but also happens to live
742 # in the parent's "natural ambient space" (this happens with
745 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, elt
)
747 natural_basis
= A
.natural_basis()
748 if elt
in natural_basis
[0].matrix_space():
749 # Thanks for nothing! Matrix spaces aren't vector
750 # spaces in Sage, so we have to figure out its
751 # natural-basis coordinates ourselves.
752 V
= VectorSpace(elt
.base_ring(), elt
.nrows()**2)
753 W
= V
.span( _mat2vec(s
) for s
in natural_basis
)
754 coords
= W
.coordinates(_mat2vec(elt
))
755 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, coords
)
757 def __pow__(self
, n
):
759 Return ``self`` raised to the power ``n``.
761 Jordan algebras are always power-associative; see for
762 example Faraut and Koranyi, Proposition II.1.2 (ii).
766 We have to override this because our superclass uses row vectors
767 instead of column vectors! We, on the other hand, assume column
772 sage: set_random_seed()
773 sage: x = random_eja().random_element()
774 sage: x.operator_matrix()*x.vector() == (x^2).vector()
777 A few examples of power-associativity::
779 sage: set_random_seed()
780 sage: x = random_eja().random_element()
781 sage: x*(x*x)*(x*x) == x^5
783 sage: (x*x)*(x*x*x) == x^5
786 We also know that powers operator-commute (Koecher, Chapter
789 sage: set_random_seed()
790 sage: x = random_eja().random_element()
791 sage: m = ZZ.random_element(0,10)
792 sage: n = ZZ.random_element(0,10)
793 sage: Lxm = (x^m).operator_matrix()
794 sage: Lxn = (x^n).operator_matrix()
795 sage: Lxm*Lxn == Lxn*Lxm
805 return A( (self
.operator_matrix()**(n
-1))*self
.vector() )
808 def apply_univariate_polynomial(self
, p
):
810 Apply the univariate polynomial ``p`` to this element.
812 A priori, SageMath won't allow us to apply a univariate
813 polynomial to an element of an EJA, because we don't know
814 that EJAs are rings (they are usually not associative). Of
815 course, we know that EJAs are power-associative, so the
816 operation is ultimately kosher. This function sidesteps
817 the CAS to get the answer we want and expect.
821 sage: R = PolynomialRing(QQ, 't')
823 sage: p = t^4 - t^3 + 5*t - 2
824 sage: J = RealCartesianProductEJA(5)
825 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
830 We should always get back an element of the algebra::
832 sage: set_random_seed()
833 sage: p = PolynomialRing(QQ, 't').random_element()
834 sage: J = random_eja()
835 sage: x = J.random_element()
836 sage: x.apply_univariate_polynomial(p) in J
840 if len(p
.variables()) > 1:
841 raise ValueError("not a univariate polynomial")
844 # Convert the coeficcients to the parent's base ring,
845 # because a priori they might live in an (unnecessarily)
846 # larger ring for which P.sum() would fail below.
847 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
848 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
851 def characteristic_polynomial(self
):
853 Return the characteristic polynomial of this element.
857 The rank of `R^3` is three, and the minimal polynomial of
858 the identity element is `(t-1)` from which it follows that
859 the characteristic polynomial should be `(t-1)^3`::
861 sage: J = RealCartesianProductEJA(3)
862 sage: J.one().characteristic_polynomial()
863 t^3 - 3*t^2 + 3*t - 1
865 Likewise, the characteristic of the zero element in the
866 rank-three algebra `R^{n}` should be `t^{3}`::
868 sage: J = RealCartesianProductEJA(3)
869 sage: J.zero().characteristic_polynomial()
872 The characteristic polynomial of an element should evaluate
873 to zero on that element::
875 sage: set_random_seed()
876 sage: x = RealCartesianProductEJA(3).random_element()
877 sage: p = x.characteristic_polynomial()
878 sage: x.apply_univariate_polynomial(p)
882 p
= self
.parent().characteristic_polynomial()
883 return p(*self
.vector())
886 def inner_product(self
, other
):
888 Return the parent algebra's inner product of myself and ``other``.
892 The inner product in the Jordan spin algebra is the usual
893 inner product on `R^n` (this example only works because the
894 basis for the Jordan algebra is the standard basis in `R^n`)::
896 sage: J = JordanSpinEJA(3)
897 sage: x = vector(QQ,[1,2,3])
898 sage: y = vector(QQ,[4,5,6])
899 sage: x.inner_product(y)
901 sage: J(x).inner_product(J(y))
904 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
905 multiplication is the usual matrix multiplication in `S^n`,
906 so the inner product of the identity matrix with itself
909 sage: J = RealSymmetricEJA(3)
910 sage: J.one().inner_product(J.one())
913 Likewise, the inner product on `C^n` is `<X,Y> =
914 Re(trace(X*Y))`, where we must necessarily take the real
915 part because the product of Hermitian matrices may not be
918 sage: J = ComplexHermitianEJA(3)
919 sage: J.one().inner_product(J.one())
922 Ditto for the quaternions::
924 sage: J = QuaternionHermitianEJA(3)
925 sage: J.one().inner_product(J.one())
930 Ensure that we can always compute an inner product, and that
931 it gives us back a real number::
933 sage: set_random_seed()
934 sage: J = random_eja()
935 sage: x = J.random_element()
936 sage: y = J.random_element()
937 sage: x.inner_product(y) in RR
943 raise TypeError("'other' must live in the same algebra")
945 return P
.inner_product(self
, other
)
948 def operator_commutes_with(self
, other
):
950 Return whether or not this element operator-commutes
955 The definition of a Jordan algebra says that any element
956 operator-commutes with its square::
958 sage: set_random_seed()
959 sage: x = random_eja().random_element()
960 sage: x.operator_commutes_with(x^2)
965 Test Lemma 1 from Chapter III of Koecher::
967 sage: set_random_seed()
968 sage: J = random_eja()
969 sage: u = J.random_element()
970 sage: v = J.random_element()
971 sage: lhs = u.operator_commutes_with(u*v)
972 sage: rhs = v.operator_commutes_with(u^2)
977 if not other
in self
.parent():
978 raise TypeError("'other' must live in the same algebra")
980 A
= self
.operator_matrix()
981 B
= other
.operator_matrix()
987 Return my determinant, the product of my eigenvalues.
991 sage: J = JordanSpinEJA(2)
992 sage: e0,e1 = J.gens()
993 sage: x = sum( J.gens() )
999 sage: J = JordanSpinEJA(3)
1000 sage: e0,e1,e2 = J.gens()
1001 sage: x = sum( J.gens() )
1007 An element is invertible if and only if its determinant is
1010 sage: set_random_seed()
1011 sage: x = random_eja().random_element()
1012 sage: x.is_invertible() == (x.det() != 0)
1018 p
= P
._charpoly
_coeff
(0)
1019 # The _charpoly_coeff function already adds the factor of
1020 # -1 to ensure that _charpoly_coeff(0) is really what
1021 # appears in front of t^{0} in the charpoly. However,
1022 # we want (-1)^r times THAT for the determinant.
1023 return ((-1)**r
)*p(*self
.vector())
1028 Return the Jordan-multiplicative inverse of this element.
1032 We appeal to the quadratic representation as in Koecher's
1033 Theorem 12 in Chapter III, Section 5.
1037 The inverse in the spin factor algebra is given in Alizadeh's
1040 sage: set_random_seed()
1041 sage: n = ZZ.random_element(1,10)
1042 sage: J = JordanSpinEJA(n)
1043 sage: x = J.random_element()
1044 sage: while x.is_zero():
1045 ....: x = J.random_element()
1046 sage: x_vec = x.vector()
1048 sage: x_bar = x_vec[1:]
1049 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
1050 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
1051 sage: x_inverse = coeff*inv_vec
1052 sage: x.inverse() == J(x_inverse)
1057 The identity element is its own inverse::
1059 sage: set_random_seed()
1060 sage: J = random_eja()
1061 sage: J.one().inverse() == J.one()
1064 If an element has an inverse, it acts like one::
1066 sage: set_random_seed()
1067 sage: J = random_eja()
1068 sage: x = J.random_element()
1069 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
1072 The inverse of the inverse is what we started with::
1074 sage: set_random_seed()
1075 sage: J = random_eja()
1076 sage: x = J.random_element()
1077 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
1080 The zero element is never invertible::
1082 sage: set_random_seed()
1083 sage: J = random_eja().zero().inverse()
1084 Traceback (most recent call last):
1086 ValueError: element is not invertible
1089 if not self
.is_invertible():
1090 raise ValueError("element is not invertible")
1092 return (~self
.quadratic_representation())(self
)
1095 def is_invertible(self
):
1097 Return whether or not this element is invertible.
1099 We can't use the superclass method because it relies on
1100 the algebra being associative.
1104 The usual way to do this is to check if the determinant is
1105 zero, but we need the characteristic polynomial for the
1106 determinant. The minimal polynomial is a lot easier to get,
1107 so we use Corollary 2 in Chapter V of Koecher to check
1108 whether or not the paren't algebra's zero element is a root
1109 of this element's minimal polynomial.
1113 The identity element is always invertible::
1115 sage: set_random_seed()
1116 sage: J = random_eja()
1117 sage: J.one().is_invertible()
1120 The zero element is never invertible::
1122 sage: set_random_seed()
1123 sage: J = random_eja()
1124 sage: J.zero().is_invertible()
1128 zero
= self
.parent().zero()
1129 p
= self
.minimal_polynomial()
1130 return not (p(zero
) == zero
)
1133 def is_nilpotent(self
):
1135 Return whether or not some power of this element is zero.
1137 The superclass method won't work unless we're in an
1138 associative algebra, and we aren't. However, we generate
1139 an assocoative subalgebra and we're nilpotent there if and
1140 only if we're nilpotent here (probably).
1144 The identity element is never nilpotent::
1146 sage: set_random_seed()
1147 sage: random_eja().one().is_nilpotent()
1150 The additive identity is always nilpotent::
1152 sage: set_random_seed()
1153 sage: random_eja().zero().is_nilpotent()
1157 # The element we're going to call "is_nilpotent()" on.
1158 # Either myself, interpreted as an element of a finite-
1159 # dimensional algebra, or an element of an associative
1163 if self
.parent().is_associative():
1164 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1166 V
= self
.span_of_powers()
1167 assoc_subalg
= self
.subalgebra_generated_by()
1168 # Mis-design warning: the basis used for span_of_powers()
1169 # and subalgebra_generated_by() must be the same, and in
1171 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1173 # Recursive call, but should work since elt lives in an
1174 # associative algebra.
1175 return elt
.is_nilpotent()
1178 def is_regular(self
):
1180 Return whether or not this is a regular element.
1184 The identity element always has degree one, but any element
1185 linearly-independent from it is regular::
1187 sage: J = JordanSpinEJA(5)
1188 sage: J.one().is_regular()
1190 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
1191 sage: for x in J.gens():
1192 ....: (J.one() + x).is_regular()
1200 return self
.degree() == self
.parent().rank()
1205 Compute the degree of this element the straightforward way
1206 according to the definition; by appending powers to a list
1207 and figuring out its dimension (that is, whether or not
1208 they're linearly dependent).
1212 sage: J = JordanSpinEJA(4)
1213 sage: J.one().degree()
1215 sage: e0,e1,e2,e3 = J.gens()
1216 sage: (e0 - e1).degree()
1219 In the spin factor algebra (of rank two), all elements that
1220 aren't multiples of the identity are regular::
1222 sage: set_random_seed()
1223 sage: n = ZZ.random_element(1,10)
1224 sage: J = JordanSpinEJA(n)
1225 sage: x = J.random_element()
1226 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
1230 return self
.span_of_powers().dimension()
1233 def left_matrix(self
):
1235 Our parent class defines ``left_matrix`` and ``matrix``
1236 methods whose names are misleading. We don't want them.
1238 raise NotImplementedError("use operator_matrix() instead")
1240 matrix
= left_matrix
1243 def minimal_polynomial(self
):
1245 Return the minimal polynomial of this element,
1246 as a function of the variable `t`.
1250 We restrict ourselves to the associative subalgebra
1251 generated by this element, and then return the minimal
1252 polynomial of this element's operator matrix (in that
1253 subalgebra). This works by Baes Proposition 2.3.16.
1257 The minimal polynomial of the identity and zero elements are
1260 sage: set_random_seed()
1261 sage: J = random_eja()
1262 sage: J.one().minimal_polynomial()
1264 sage: J.zero().minimal_polynomial()
1267 The degree of an element is (by one definition) the degree
1268 of its minimal polynomial::
1270 sage: set_random_seed()
1271 sage: x = random_eja().random_element()
1272 sage: x.degree() == x.minimal_polynomial().degree()
1275 The minimal polynomial and the characteristic polynomial coincide
1276 and are known (see Alizadeh, Example 11.11) for all elements of
1277 the spin factor algebra that aren't scalar multiples of the
1280 sage: set_random_seed()
1281 sage: n = ZZ.random_element(2,10)
1282 sage: J = JordanSpinEJA(n)
1283 sage: y = J.random_element()
1284 sage: while y == y.coefficient(0)*J.one():
1285 ....: y = J.random_element()
1286 sage: y0 = y.vector()[0]
1287 sage: y_bar = y.vector()[1:]
1288 sage: actual = y.minimal_polynomial()
1289 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1290 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1291 sage: bool(actual == expected)
1294 The minimal polynomial should always kill its element::
1296 sage: set_random_seed()
1297 sage: x = random_eja().random_element()
1298 sage: p = x.minimal_polynomial()
1299 sage: x.apply_univariate_polynomial(p)
1303 V
= self
.span_of_powers()
1304 assoc_subalg
= self
.subalgebra_generated_by()
1305 # Mis-design warning: the basis used for span_of_powers()
1306 # and subalgebra_generated_by() must be the same, and in
1308 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1310 # We get back a symbolic polynomial in 'x' but want a real
1311 # polynomial in 't'.
1312 p_of_x
= elt
.operator_matrix().minimal_polynomial()
1313 return p_of_x
.change_variable_name('t')
1316 def natural_representation(self
):
1318 Return a more-natural representation of this element.
1320 Every finite-dimensional Euclidean Jordan Algebra is a
1321 direct sum of five simple algebras, four of which comprise
1322 Hermitian matrices. This method returns the original
1323 "natural" representation of this element as a Hermitian
1324 matrix, if it has one. If not, you get the usual representation.
1328 sage: J = ComplexHermitianEJA(3)
1331 sage: J.one().natural_representation()
1341 sage: J = QuaternionHermitianEJA(3)
1344 sage: J.one().natural_representation()
1345 [1 0 0 0 0 0 0 0 0 0 0 0]
1346 [0 1 0 0 0 0 0 0 0 0 0 0]
1347 [0 0 1 0 0 0 0 0 0 0 0 0]
1348 [0 0 0 1 0 0 0 0 0 0 0 0]
1349 [0 0 0 0 1 0 0 0 0 0 0 0]
1350 [0 0 0 0 0 1 0 0 0 0 0 0]
1351 [0 0 0 0 0 0 1 0 0 0 0 0]
1352 [0 0 0 0 0 0 0 1 0 0 0 0]
1353 [0 0 0 0 0 0 0 0 1 0 0 0]
1354 [0 0 0 0 0 0 0 0 0 1 0 0]
1355 [0 0 0 0 0 0 0 0 0 0 1 0]
1356 [0 0 0 0 0 0 0 0 0 0 0 1]
1359 B
= self
.parent().natural_basis()
1360 W
= B
[0].matrix_space()
1361 return W
.linear_combination(zip(self
.vector(), B
))
1366 Return the left-multiplication-by-this-element
1367 operator on the ambient algebra.
1371 sage: set_random_seed()
1372 sage: J = random_eja()
1373 sage: x = J.random_element()
1374 sage: y = J.random_element()
1375 sage: x.operator()(y) == x*y
1377 sage: y.operator()(x) == x*y
1382 return FiniteDimensionalEuclideanJordanAlgebraMorphism(
1384 self
.operator_matrix() )
1388 def operator_matrix(self
):
1390 Return the matrix that represents left- (or right-)
1391 multiplication by this element in the parent algebra.
1393 We implement this ourselves to work around the fact that
1394 our parent class represents everything with row vectors.
1398 Test the first polarization identity from my notes, Koecher Chapter
1399 III, or from Baes (2.3)::
1401 sage: set_random_seed()
1402 sage: J = random_eja()
1403 sage: x = J.random_element()
1404 sage: y = J.random_element()
1405 sage: Lx = x.operator_matrix()
1406 sage: Ly = y.operator_matrix()
1407 sage: Lxx = (x*x).operator_matrix()
1408 sage: Lxy = (x*y).operator_matrix()
1409 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
1412 Test the second polarization identity from my notes or from
1415 sage: set_random_seed()
1416 sage: J = random_eja()
1417 sage: x = J.random_element()
1418 sage: y = J.random_element()
1419 sage: z = J.random_element()
1420 sage: Lx = x.operator_matrix()
1421 sage: Ly = y.operator_matrix()
1422 sage: Lz = z.operator_matrix()
1423 sage: Lzy = (z*y).operator_matrix()
1424 sage: Lxy = (x*y).operator_matrix()
1425 sage: Lxz = (x*z).operator_matrix()
1426 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
1429 Test the third polarization identity from my notes or from
1432 sage: set_random_seed()
1433 sage: J = random_eja()
1434 sage: u = J.random_element()
1435 sage: y = J.random_element()
1436 sage: z = J.random_element()
1437 sage: Lu = u.operator_matrix()
1438 sage: Ly = y.operator_matrix()
1439 sage: Lz = z.operator_matrix()
1440 sage: Lzy = (z*y).operator_matrix()
1441 sage: Luy = (u*y).operator_matrix()
1442 sage: Luz = (u*z).operator_matrix()
1443 sage: Luyz = (u*(y*z)).operator_matrix()
1444 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
1445 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
1446 sage: bool(lhs == rhs)
1450 fda_elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1451 return fda_elt
.matrix().transpose()
1454 def quadratic_representation(self
, other
=None):
1456 Return the quadratic representation of this element.
1460 The explicit form in the spin factor algebra is given by
1461 Alizadeh's Example 11.12::
1463 sage: set_random_seed()
1464 sage: n = ZZ.random_element(1,10)
1465 sage: J = JordanSpinEJA(n)
1466 sage: x = J.random_element()
1467 sage: x_vec = x.vector()
1469 sage: x_bar = x_vec[1:]
1470 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1471 sage: B = 2*x0*x_bar.row()
1472 sage: C = 2*x0*x_bar.column()
1473 sage: D = identity_matrix(QQ, n-1)
1474 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1475 sage: D = D + 2*x_bar.tensor_product(x_bar)
1476 sage: Q = block_matrix(2,2,[A,B,C,D])
1477 sage: Q == x.quadratic_representation().matrix()
1480 Test all of the properties from Theorem 11.2 in Alizadeh::
1482 sage: set_random_seed()
1483 sage: J = random_eja()
1484 sage: x = J.random_element()
1485 sage: y = J.random_element()
1486 sage: Lx = x.operator()
1487 sage: Lxx = (x*x).operator()
1488 sage: Qx = x.quadratic_representation()
1489 sage: Qy = y.quadratic_representation()
1490 sage: Qxy = x.quadratic_representation(y)
1491 sage: Qex = J.one().quadratic_representation(x)
1492 sage: n = ZZ.random_element(10)
1493 sage: Qxn = (x^n).quadratic_representation()
1497 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1502 sage: alpha = QQ.random_element()
1503 sage: (alpha*x).quadratic_representation() == (alpha^2)*Qx
1508 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1511 sage: not x.is_invertible() or (
1514 ....: x.inverse().quadratic_representation() )
1517 sage: Qxy(J.one()) == x*y
1522 sage: not x.is_invertible() or (
1523 ....: x.quadratic_representation(x.inverse())*Qx
1524 ....: == Qx*x.quadratic_representation(x.inverse()) )
1527 sage: not x.is_invertible() or (
1528 ....: x.quadratic_representation(x.inverse())*Qx
1530 ....: 2*x.operator()*Qex - Qx )
1533 sage: 2*x.operator()*Qex - Qx == Lxx
1538 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1548 sage: not x.is_invertible() or (
1549 ....: Qx*x.inverse().operator() == Lx )
1554 sage: not x.operator_commutes_with(y) or (
1555 ....: Qx(y)^n == Qxn(y^n) )
1561 elif not other
in self
.parent():
1562 raise TypeError("'other' must live in the same algebra")
1565 M
= other
.operator()
1566 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1569 def span_of_powers(self
):
1571 Return the vector space spanned by successive powers of
1574 # The dimension of the subalgebra can't be greater than
1575 # the big algebra, so just put everything into a list
1576 # and let span() get rid of the excess.
1578 # We do the extra ambient_vector_space() in case we're messing
1579 # with polynomials and the direct parent is a module.
1580 V
= self
.vector().parent().ambient_vector_space()
1581 return V
.span( (self
**d
).vector() for d
in xrange(V
.dimension()) )
1584 def subalgebra_generated_by(self
):
1586 Return the associative subalgebra of the parent EJA generated
1591 sage: set_random_seed()
1592 sage: x = random_eja().random_element()
1593 sage: x.subalgebra_generated_by().is_associative()
1596 Squaring in the subalgebra should be the same thing as
1597 squaring in the superalgebra::
1599 sage: set_random_seed()
1600 sage: x = random_eja().random_element()
1601 sage: u = x.subalgebra_generated_by().random_element()
1602 sage: u.operator_matrix()*u.vector() == (u**2).vector()
1606 # First get the subspace spanned by the powers of myself...
1607 V
= self
.span_of_powers()
1608 F
= self
.base_ring()
1610 # Now figure out the entries of the right-multiplication
1611 # matrix for the successive basis elements b0, b1,... of
1614 for b_right
in V
.basis():
1615 eja_b_right
= self
.parent()(b_right
)
1617 # The first row of the right-multiplication matrix by
1618 # b1 is what we get if we apply that matrix to b1. The
1619 # second row of the right multiplication matrix by b1
1620 # is what we get when we apply that matrix to b2...
1622 # IMPORTANT: this assumes that all vectors are COLUMN
1623 # vectors, unlike our superclass (which uses row vectors).
1624 for b_left
in V
.basis():
1625 eja_b_left
= self
.parent()(b_left
)
1626 # Multiply in the original EJA, but then get the
1627 # coordinates from the subalgebra in terms of its
1629 this_row
= V
.coordinates((eja_b_left
*eja_b_right
).vector())
1630 b_right_rows
.append(this_row
)
1631 b_right_matrix
= matrix(F
, b_right_rows
)
1632 mats
.append(b_right_matrix
)
1634 # It's an algebra of polynomials in one element, and EJAs
1635 # are power-associative.
1637 # TODO: choose generator names intelligently.
1638 return FiniteDimensionalEuclideanJordanAlgebra(F
, mats
, assume_associative
=True, names
='f')
1641 def subalgebra_idempotent(self
):
1643 Find an idempotent in the associative subalgebra I generate
1644 using Proposition 2.3.5 in Baes.
1648 sage: set_random_seed()
1649 sage: J = random_eja()
1650 sage: x = J.random_element()
1651 sage: while x.is_nilpotent():
1652 ....: x = J.random_element()
1653 sage: c = x.subalgebra_idempotent()
1658 if self
.is_nilpotent():
1659 raise ValueError("this only works with non-nilpotent elements!")
1661 V
= self
.span_of_powers()
1662 J
= self
.subalgebra_generated_by()
1663 # Mis-design warning: the basis used for span_of_powers()
1664 # and subalgebra_generated_by() must be the same, and in
1666 u
= J(V
.coordinates(self
.vector()))
1668 # The image of the matrix of left-u^m-multiplication
1669 # will be minimal for some natural number s...
1671 minimal_dim
= V
.dimension()
1672 for i
in xrange(1, V
.dimension()):
1673 this_dim
= (u
**i
).operator_matrix().image().dimension()
1674 if this_dim
< minimal_dim
:
1675 minimal_dim
= this_dim
1678 # Now minimal_matrix should correspond to the smallest
1679 # non-zero subspace in Baes's (or really, Koecher's)
1682 # However, we need to restrict the matrix to work on the
1683 # subspace... or do we? Can't we just solve, knowing that
1684 # A(c) = u^(s+1) should have a solution in the big space,
1687 # Beware, solve_right() means that we're using COLUMN vectors.
1688 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1690 A
= u_next
.operator_matrix()
1691 c_coordinates
= A
.solve_right(u_next
.vector())
1693 # Now c_coordinates is the idempotent we want, but it's in
1694 # the coordinate system of the subalgebra.
1696 # We need the basis for J, but as elements of the parent algebra.
1698 basis
= [self
.parent(v
) for v
in V
.basis()]
1699 return self
.parent().linear_combination(zip(c_coordinates
, basis
))
1704 Return my trace, the sum of my eigenvalues.
1708 sage: J = JordanSpinEJA(3)
1709 sage: x = sum(J.gens())
1715 sage: J = RealCartesianProductEJA(5)
1716 sage: J.one().trace()
1721 The trace of an element is a real number::
1723 sage: set_random_seed()
1724 sage: J = random_eja()
1725 sage: J.random_element().trace() in J.base_ring()
1731 p
= P
._charpoly
_coeff
(r
-1)
1732 # The _charpoly_coeff function already adds the factor of
1733 # -1 to ensure that _charpoly_coeff(r-1) is really what
1734 # appears in front of t^{r-1} in the charpoly. However,
1735 # we want the negative of THAT for the trace.
1736 return -p(*self
.vector())
1739 def trace_inner_product(self
, other
):
1741 Return the trace inner product of myself and ``other``.
1745 The trace inner product is commutative::
1747 sage: set_random_seed()
1748 sage: J = random_eja()
1749 sage: x = J.random_element(); y = J.random_element()
1750 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1753 The trace inner product is bilinear::
1755 sage: set_random_seed()
1756 sage: J = random_eja()
1757 sage: x = J.random_element()
1758 sage: y = J.random_element()
1759 sage: z = J.random_element()
1760 sage: a = QQ.random_element();
1761 sage: actual = (a*(x+z)).trace_inner_product(y)
1762 sage: expected = ( a*x.trace_inner_product(y) +
1763 ....: a*z.trace_inner_product(y) )
1764 sage: actual == expected
1766 sage: actual = x.trace_inner_product(a*(y+z))
1767 sage: expected = ( a*x.trace_inner_product(y) +
1768 ....: a*x.trace_inner_product(z) )
1769 sage: actual == expected
1772 The trace inner product satisfies the compatibility
1773 condition in the definition of a Euclidean Jordan algebra::
1775 sage: set_random_seed()
1776 sage: J = random_eja()
1777 sage: x = J.random_element()
1778 sage: y = J.random_element()
1779 sage: z = J.random_element()
1780 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1784 if not other
in self
.parent():
1785 raise TypeError("'other' must live in the same algebra")
1787 return (self
*other
).trace()
1790 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1792 Return the Euclidean Jordan Algebra corresponding to the set
1793 `R^n` under the Hadamard product.
1795 Note: this is nothing more than the Cartesian product of ``n``
1796 copies of the spin algebra. Once Cartesian product algebras
1797 are implemented, this can go.
1801 This multiplication table can be verified by hand::
1803 sage: J = RealCartesianProductEJA(3)
1804 sage: e0,e1,e2 = J.gens()
1820 def __classcall_private__(cls
, n
, field
=QQ
):
1821 # The FiniteDimensionalAlgebra constructor takes a list of
1822 # matrices, the ith representing right multiplication by the ith
1823 # basis element in the vector space. So if e_1 = (1,0,0), then
1824 # right (Hadamard) multiplication of x by e_1 picks out the first
1825 # component of x; and likewise for the ith basis element e_i.
1826 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
1827 for i
in xrange(n
) ]
1829 fdeja
= super(RealCartesianProductEJA
, cls
)
1830 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=n
)
1832 def inner_product(self
, x
, y
):
1833 return _usual_ip(x
,y
)
1838 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1842 For now, we choose a random natural number ``n`` (greater than zero)
1843 and then give you back one of the following:
1845 * The cartesian product of the rational numbers ``n`` times; this is
1846 ``QQ^n`` with the Hadamard product.
1848 * The Jordan spin algebra on ``QQ^n``.
1850 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1853 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1854 in the space of ``2n``-by-``2n`` real symmetric matrices.
1856 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1857 in the space of ``4n``-by-``4n`` real symmetric matrices.
1859 Later this might be extended to return Cartesian products of the
1865 Euclidean Jordan algebra of degree...
1869 # The max_n component lets us choose different upper bounds on the
1870 # value "n" that gets passed to the constructor. This is needed
1871 # because e.g. R^{10} is reasonable to test, while the Hermitian
1872 # 10-by-10 quaternion matrices are not.
1873 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
1875 (RealSymmetricEJA
, 5),
1876 (ComplexHermitianEJA
, 4),
1877 (QuaternionHermitianEJA
, 3)])
1878 n
= ZZ
.random_element(1, max_n
)
1879 return constructor(n
, field
=QQ
)
1883 def _real_symmetric_basis(n
, field
=QQ
):
1885 Return a basis for the space of real symmetric n-by-n matrices.
1887 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1891 for j
in xrange(i
+1):
1892 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1896 # Beware, orthogonal but not normalized!
1897 Sij
= Eij
+ Eij
.transpose()
1902 def _complex_hermitian_basis(n
, field
=QQ
):
1904 Returns a basis for the space of complex Hermitian n-by-n matrices.
1908 sage: set_random_seed()
1909 sage: n = ZZ.random_element(1,5)
1910 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1914 F
= QuadraticField(-1, 'I')
1917 # This is like the symmetric case, but we need to be careful:
1919 # * We want conjugate-symmetry, not just symmetry.
1920 # * The diagonal will (as a result) be real.
1924 for j
in xrange(i
+1):
1925 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1927 Sij
= _embed_complex_matrix(Eij
)
1930 # Beware, orthogonal but not normalized! The second one
1931 # has a minus because it's conjugated.
1932 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
1934 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
1939 def _quaternion_hermitian_basis(n
, field
=QQ
):
1941 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1945 sage: set_random_seed()
1946 sage: n = ZZ.random_element(1,5)
1947 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1951 Q
= QuaternionAlgebra(QQ
,-1,-1)
1954 # This is like the symmetric case, but we need to be careful:
1956 # * We want conjugate-symmetry, not just symmetry.
1957 # * The diagonal will (as a result) be real.
1961 for j
in xrange(i
+1):
1962 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1964 Sij
= _embed_quaternion_matrix(Eij
)
1967 # Beware, orthogonal but not normalized! The second,
1968 # third, and fourth ones have a minus because they're
1970 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
1972 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
1974 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
1976 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
1982 return vector(m
.base_ring(), m
.list())
1985 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
1987 def _multiplication_table_from_matrix_basis(basis
):
1989 At least three of the five simple Euclidean Jordan algebras have the
1990 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1991 multiplication on the right is matrix multiplication. Given a basis
1992 for the underlying matrix space, this function returns a
1993 multiplication table (obtained by looping through the basis
1994 elements) for an algebra of those matrices. A reordered copy
1995 of the basis is also returned to work around the fact that
1996 the ``span()`` in this function will change the order of the basis
1997 from what we think it is, to... something else.
1999 # In S^2, for example, we nominally have four coordinates even
2000 # though the space is of dimension three only. The vector space V
2001 # is supposed to hold the entire long vector, and the subspace W
2002 # of V will be spanned by the vectors that arise from symmetric
2003 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
2004 field
= basis
[0].base_ring()
2005 dimension
= basis
[0].nrows()
2007 V
= VectorSpace(field
, dimension
**2)
2008 W
= V
.span( _mat2vec(s
) for s
in basis
)
2010 # Taking the span above reorders our basis (thanks, jerk!) so we
2011 # need to put our "matrix basis" in the same order as the
2012 # (reordered) vector basis.
2013 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
2017 # Brute force the multiplication-by-s matrix by looping
2018 # through all elements of the basis and doing the computation
2019 # to find out what the corresponding row should be. BEWARE:
2020 # these multiplication tables won't be symmetric! It therefore
2021 # becomes REALLY IMPORTANT that the underlying algebra
2022 # constructor uses ROW vectors and not COLUMN vectors. That's
2023 # why we're computing rows here and not columns.
2026 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
2027 Q_rows
.append(W
.coordinates(this_row
))
2028 Q
= matrix(field
, W
.dimension(), Q_rows
)
2034 def _embed_complex_matrix(M
):
2036 Embed the n-by-n complex matrix ``M`` into the space of real
2037 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2038 bi` to the block matrix ``[[a,b],[-b,a]]``.
2042 sage: F = QuadraticField(-1,'i')
2043 sage: x1 = F(4 - 2*i)
2044 sage: x2 = F(1 + 2*i)
2047 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2048 sage: _embed_complex_matrix(M)
2057 Embedding is a homomorphism (isomorphism, in fact)::
2059 sage: set_random_seed()
2060 sage: n = ZZ.random_element(5)
2061 sage: F = QuadraticField(-1, 'i')
2062 sage: X = random_matrix(F, n)
2063 sage: Y = random_matrix(F, n)
2064 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
2065 sage: expected = _embed_complex_matrix(X*Y)
2066 sage: actual == expected
2072 raise ValueError("the matrix 'M' must be square")
2073 field
= M
.base_ring()
2078 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
2080 # We can drop the imaginaries here.
2081 return block_matrix(field
.base_ring(), n
, blocks
)
2084 def _unembed_complex_matrix(M
):
2086 The inverse of _embed_complex_matrix().
2090 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2091 ....: [-2, 1, -4, 3],
2092 ....: [ 9, 10, 11, 12],
2093 ....: [-10, 9, -12, 11] ])
2094 sage: _unembed_complex_matrix(A)
2096 [ 10*i + 9 12*i + 11]
2100 Unembedding is the inverse of embedding::
2102 sage: set_random_seed()
2103 sage: F = QuadraticField(-1, 'i')
2104 sage: M = random_matrix(F, 3)
2105 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
2111 raise ValueError("the matrix 'M' must be square")
2112 if not n
.mod(2).is_zero():
2113 raise ValueError("the matrix 'M' must be a complex embedding")
2115 F
= QuadraticField(-1, 'i')
2118 # Go top-left to bottom-right (reading order), converting every
2119 # 2-by-2 block we see to a single complex element.
2121 for k
in xrange(n
/2):
2122 for j
in xrange(n
/2):
2123 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
2124 if submat
[0,0] != submat
[1,1]:
2125 raise ValueError('bad on-diagonal submatrix')
2126 if submat
[0,1] != -submat
[1,0]:
2127 raise ValueError('bad off-diagonal submatrix')
2128 z
= submat
[0,0] + submat
[0,1]*i
2131 return matrix(F
, n
/2, elements
)
2134 def _embed_quaternion_matrix(M
):
2136 Embed the n-by-n quaternion matrix ``M`` into the space of real
2137 matrices of size 4n-by-4n by first sending each quaternion entry
2138 `z = a + bi + cj + dk` to the block-complex matrix
2139 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
2144 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2145 sage: i,j,k = Q.gens()
2146 sage: x = 1 + 2*i + 3*j + 4*k
2147 sage: M = matrix(Q, 1, [[x]])
2148 sage: _embed_quaternion_matrix(M)
2154 Embedding is a homomorphism (isomorphism, in fact)::
2156 sage: set_random_seed()
2157 sage: n = ZZ.random_element(5)
2158 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2159 sage: X = random_matrix(Q, n)
2160 sage: Y = random_matrix(Q, n)
2161 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
2162 sage: expected = _embed_quaternion_matrix(X*Y)
2163 sage: actual == expected
2167 quaternions
= M
.base_ring()
2170 raise ValueError("the matrix 'M' must be square")
2172 F
= QuadraticField(-1, 'i')
2177 t
= z
.coefficient_tuple()
2182 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2183 [-c
+ d
*i
, a
- b
*i
]])
2184 blocks
.append(_embed_complex_matrix(cplx_matrix
))
2186 # We should have real entries by now, so use the realest field
2187 # we've got for the return value.
2188 return block_matrix(quaternions
.base_ring(), n
, blocks
)
2191 def _unembed_quaternion_matrix(M
):
2193 The inverse of _embed_quaternion_matrix().
2197 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2198 ....: [-2, 1, -4, 3],
2199 ....: [-3, 4, 1, -2],
2200 ....: [-4, -3, 2, 1]])
2201 sage: _unembed_quaternion_matrix(M)
2202 [1 + 2*i + 3*j + 4*k]
2206 Unembedding is the inverse of embedding::
2208 sage: set_random_seed()
2209 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2210 sage: M = random_matrix(Q, 3)
2211 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
2217 raise ValueError("the matrix 'M' must be square")
2218 if not n
.mod(4).is_zero():
2219 raise ValueError("the matrix 'M' must be a complex embedding")
2221 Q
= QuaternionAlgebra(QQ
,-1,-1)
2224 # Go top-left to bottom-right (reading order), converting every
2225 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2228 for l
in xrange(n
/4):
2229 for m
in xrange(n
/4):
2230 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
2231 if submat
[0,0] != submat
[1,1].conjugate():
2232 raise ValueError('bad on-diagonal submatrix')
2233 if submat
[0,1] != -submat
[1,0].conjugate():
2234 raise ValueError('bad off-diagonal submatrix')
2235 z
= submat
[0,0].real() + submat
[0,0].imag()*i
2236 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
2239 return matrix(Q
, n
/4, elements
)
2242 # The usual inner product on R^n.
2244 return x
.vector().inner_product(y
.vector())
2246 # The inner product used for the real symmetric simple EJA.
2247 # We keep it as a separate function because e.g. the complex
2248 # algebra uses the same inner product, except divided by 2.
2249 def _matrix_ip(X
,Y
):
2250 X_mat
= X
.natural_representation()
2251 Y_mat
= Y
.natural_representation()
2252 return (X_mat
*Y_mat
).trace()
2255 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2257 The rank-n simple EJA consisting of real symmetric n-by-n
2258 matrices, the usual symmetric Jordan product, and the trace inner
2259 product. It has dimension `(n^2 + n)/2` over the reals.
2263 sage: J = RealSymmetricEJA(2)
2264 sage: e0, e1, e2 = J.gens()
2274 The degree of this algebra is `(n^2 + n) / 2`::
2276 sage: set_random_seed()
2277 sage: n = ZZ.random_element(1,5)
2278 sage: J = RealSymmetricEJA(n)
2279 sage: J.degree() == (n^2 + n)/2
2282 The Jordan multiplication is what we think it is::
2284 sage: set_random_seed()
2285 sage: n = ZZ.random_element(1,5)
2286 sage: J = RealSymmetricEJA(n)
2287 sage: x = J.random_element()
2288 sage: y = J.random_element()
2289 sage: actual = (x*y).natural_representation()
2290 sage: X = x.natural_representation()
2291 sage: Y = y.natural_representation()
2292 sage: expected = (X*Y + Y*X)/2
2293 sage: actual == expected
2295 sage: J(expected) == x*y
2300 def __classcall_private__(cls
, n
, field
=QQ
):
2301 S
= _real_symmetric_basis(n
, field
=field
)
2302 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2304 fdeja
= super(RealSymmetricEJA
, cls
)
2305 return fdeja
.__classcall
_private
__(cls
,
2311 def inner_product(self
, x
, y
):
2312 return _matrix_ip(x
,y
)
2315 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2317 The rank-n simple EJA consisting of complex Hermitian n-by-n
2318 matrices over the real numbers, the usual symmetric Jordan product,
2319 and the real-part-of-trace inner product. It has dimension `n^2` over
2324 The degree of this algebra is `n^2`::
2326 sage: set_random_seed()
2327 sage: n = ZZ.random_element(1,5)
2328 sage: J = ComplexHermitianEJA(n)
2329 sage: J.degree() == n^2
2332 The Jordan multiplication is what we think it is::
2334 sage: set_random_seed()
2335 sage: n = ZZ.random_element(1,5)
2336 sage: J = ComplexHermitianEJA(n)
2337 sage: x = J.random_element()
2338 sage: y = J.random_element()
2339 sage: actual = (x*y).natural_representation()
2340 sage: X = x.natural_representation()
2341 sage: Y = y.natural_representation()
2342 sage: expected = (X*Y + Y*X)/2
2343 sage: actual == expected
2345 sage: J(expected) == x*y
2350 def __classcall_private__(cls
, n
, field
=QQ
):
2351 S
= _complex_hermitian_basis(n
)
2352 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2354 fdeja
= super(ComplexHermitianEJA
, cls
)
2355 return fdeja
.__classcall
_private
__(cls
,
2361 def inner_product(self
, x
, y
):
2362 # Since a+bi on the diagonal is represented as
2367 # we'll double-count the "a" entries if we take the trace of
2369 return _matrix_ip(x
,y
)/2
2372 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2374 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2375 matrices, the usual symmetric Jordan product, and the
2376 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2381 The degree of this algebra is `n^2`::
2383 sage: set_random_seed()
2384 sage: n = ZZ.random_element(1,5)
2385 sage: J = QuaternionHermitianEJA(n)
2386 sage: J.degree() == 2*(n^2) - n
2389 The Jordan multiplication is what we think it is::
2391 sage: set_random_seed()
2392 sage: n = ZZ.random_element(1,5)
2393 sage: J = QuaternionHermitianEJA(n)
2394 sage: x = J.random_element()
2395 sage: y = J.random_element()
2396 sage: actual = (x*y).natural_representation()
2397 sage: X = x.natural_representation()
2398 sage: Y = y.natural_representation()
2399 sage: expected = (X*Y + Y*X)/2
2400 sage: actual == expected
2402 sage: J(expected) == x*y
2407 def __classcall_private__(cls
, n
, field
=QQ
):
2408 S
= _quaternion_hermitian_basis(n
)
2409 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2411 fdeja
= super(QuaternionHermitianEJA
, cls
)
2412 return fdeja
.__classcall
_private
__(cls
,
2418 def inner_product(self
, x
, y
):
2419 # Since a+bi+cj+dk on the diagonal is represented as
2421 # a + bi +cj + dk = [ a b c d]
2426 # we'll quadruple-count the "a" entries if we take the trace of
2428 return _matrix_ip(x
,y
)/4
2431 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2433 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2434 with the usual inner product and jordan product ``x*y =
2435 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2440 This multiplication table can be verified by hand::
2442 sage: J = JordanSpinEJA(4)
2443 sage: e0,e1,e2,e3 = J.gens()
2461 def __classcall_private__(cls
, n
, field
=QQ
):
2463 id_matrix
= identity_matrix(field
, n
)
2465 ei
= id_matrix
.column(i
)
2466 Qi
= zero_matrix(field
, n
)
2468 Qi
.set_column(0, ei
)
2469 Qi
+= diagonal_matrix(n
, [ei
[0]]*n
)
2470 # The addition of the diagonal matrix adds an extra ei[0] in the
2471 # upper-left corner of the matrix.
2472 Qi
[0,0] = Qi
[0,0] * ~
field(2)
2475 # The rank of the spin algebra is two, unless we're in a
2476 # one-dimensional ambient space (because the rank is bounded by
2477 # the ambient dimension).
2478 fdeja
= super(JordanSpinEJA
, cls
)
2479 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=min(n
,2))
2481 def inner_product(self
, x
, y
):
2482 return _usual_ip(x
,y
)