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
14 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra_morphism
import FiniteDimensionalAlgebraMorphism
17 class FiniteDimensionalEuclideanJordanAlgebraMorphism(FiniteDimensionalAlgebraMorphism
):
19 A linear map between two finite-dimensional EJAs.
21 This is a very thin wrapper around FiniteDimensionalAlgebraMorphism
22 that does only a few things:
24 1. Avoids the ``unitary`` and ``check`` arguments to the constructor
25 that will always be ``False``. This is necessary because these
26 are homomorphisms with respect to ADDITION, but the SageMath
27 machinery wants to check that they're homomorphisms with respect
28 to (Jordan) MULTIPLICATION. That obviously doesn't work.
30 2. Inputs and outputs the underlying matrix with respect to COLUMN
31 vectors, unlike the parent class.
33 3. Allows us to add, multiply (compose), and invert morphisms in
36 If this seems a bit heavyweight, it is. I would have been happy to
37 use a the ring morphism that underlies the finite-dimensional
38 algebra morphism, but they don't seem to be callable on elements of
39 our EJA, and you can't add/multiply/invert them.
42 def __add__(self
, other
):
44 Add two EJA morphisms in the obvious way.
48 sage: J = RealSymmetricEJA(3)
51 sage: x.operator() + y.operator()
52 Morphism from Euclidean Jordan algebra of degree 6 over Rational
53 Field to Euclidean Jordan algebra of degree 6 over Rational Field
64 sage: set_random_seed()
65 sage: J = random_eja()
66 sage: x = J.random_element()
67 sage: y = J.random_element()
68 sage: (x.operator() + y.operator()) in J.Hom(J)
74 raise ValueError("summands must live in the same space")
76 return FiniteDimensionalEuclideanJordanAlgebraMorphism(
78 self
.matrix() + other
.matrix() )
81 def __init__(self
, parent
, f
):
82 FiniteDimensionalAlgebraMorphism
.__init
__(self
,
93 sage: J = RealSymmetricEJA(2)
94 sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
95 sage: x.is_invertible()
98 Morphism from Euclidean Jordan algebra of degree 3 over Rational
99 Field to Euclidean Jordan algebra of degree 3 over Rational Field
104 sage: x.operator_matrix().inverse()
111 sage: set_random_seed()
112 sage: J = random_eja()
113 sage: x = J.random_element()
114 sage: not x.is_invertible() or (
115 ....: (~x.operator()).matrix() == x.operator_matrix().inverse() )
120 if not A
.is_invertible():
121 raise ValueError("morphism is not invertible")
124 return FiniteDimensionalEuclideanJordanAlgebraMorphism(self
.parent(),
127 def __mul__(self
, other
):
129 Compose two EJA morphisms using multiplicative notation.
133 sage: J = RealSymmetricEJA(3)
136 sage: x.operator() * y.operator()
137 Morphism from Euclidean Jordan algebra of degree 6 over Rational
138 Field to Euclidean Jordan algebra of degree 6 over Rational Field
149 sage: set_random_seed()
150 sage: J = random_eja()
151 sage: x = J.random_element()
152 sage: y = J.random_element()
153 sage: (x.operator() * y.operator()) in J.Hom(J)
157 if not other
.codomain() is self
.domain():
158 raise ValueError("(co)domains must agree for composition")
160 return FiniteDimensionalEuclideanJordanAlgebraMorphism(
162 self
.matrix()*other
.matrix() )
167 We override only the representation that is shown to the user,
168 because we want the matrix to be with respect to COLUMN vectors.
172 Ensure that we see the transpose of the underlying matrix object:
174 sage: J = RealSymmetricEJA(3)
175 sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
176 sage: L = x.operator()
178 Morphism from Euclidean Jordan algebra of degree 6 over Rational
179 Field to Euclidean Jordan algebra of degree 6 over Rational Field
196 return "Morphism from {} to {} given by matrix\n{}".format(
197 self
.domain(), self
.codomain(), self
.matrix())
201 Return the matrix of this morphism with respect to a left-action
204 return FiniteDimensionalAlgebraMorphism
.matrix(self
).transpose()
207 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra
):
209 def __classcall_private__(cls
,
213 assume_associative
=False,
218 mult_table
= [b
.base_extend(field
) for b
in mult_table
]
221 if not (is_Matrix(b
) and b
.dimensions() == (n
, n
)):
222 raise ValueError("input is not a multiplication table")
223 mult_table
= tuple(mult_table
)
225 cat
= MagmaticAlgebras(field
).FiniteDimensional().WithBasis()
226 cat
.or_subcategory(category
)
227 if assume_associative
:
228 cat
= cat
.Associative()
230 names
= normalize_names(n
, names
)
232 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, cls
)
233 return fda
.__classcall
__(cls
,
236 assume_associative
=assume_associative
,
240 natural_basis
=natural_basis
)
247 assume_associative
=False,
254 By definition, Jordan multiplication commutes::
256 sage: set_random_seed()
257 sage: J = random_eja()
258 sage: x = J.random_element()
259 sage: y = J.random_element()
265 self
._natural
_basis
= natural_basis
266 self
._multiplication
_table
= mult_table
267 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
276 Return a string representation of ``self``.
278 fmt
= "Euclidean Jordan algebra of degree {} over {}"
279 return fmt
.format(self
.degree(), self
.base_ring())
282 def _a_regular_element(self
):
284 Guess a regular element. Needed to compute the basis for our
285 characteristic polynomial coefficients.
288 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
289 if not z
.is_regular():
290 raise ValueError("don't know a regular element")
295 def _charpoly_basis_space(self
):
297 Return the vector space spanned by the basis used in our
298 characteristic polynomial coefficients. This is used not only to
299 compute those coefficients, but also any time we need to
300 evaluate the coefficients (like when we compute the trace or
303 z
= self
._a
_regular
_element
()
304 V
= z
.vector().parent().ambient_vector_space()
305 V1
= V
.span_of_basis( (z
**k
).vector() for k
in range(self
.rank()) )
306 b
= (V1
.basis() + V1
.complement().basis())
307 return V
.span_of_basis(b
)
311 def _charpoly_coeff(self
, i
):
313 Return the coefficient polynomial "a_{i}" of this algebra's
314 general characteristic polynomial.
316 Having this be a separate cached method lets us compute and
317 store the trace/determinant (a_{r-1} and a_{0} respectively)
318 separate from the entire characteristic polynomial.
320 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
321 R
= A_of_x
.base_ring()
323 # Guaranteed by theory
326 # Danger: the in-place modification is done for performance
327 # reasons (reconstructing a matrix with huge polynomial
328 # entries is slow), but I don't know how cached_method works,
329 # so it's highly possible that we're modifying some global
330 # list variable by reference, here. In other words, you
331 # probably shouldn't call this method twice on the same
332 # algebra, at the same time, in two threads
333 Ai_orig
= A_of_x
.column(i
)
334 A_of_x
.set_column(i
,xr
)
335 numerator
= A_of_x
.det()
336 A_of_x
.set_column(i
,Ai_orig
)
338 # We're relying on the theory here to ensure that each a_i is
339 # indeed back in R, and the added negative signs are to make
340 # the whole charpoly expression sum to zero.
341 return R(-numerator
/detA
)
345 def _charpoly_matrix_system(self
):
347 Compute the matrix whose entries A_ij are polynomials in
348 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
349 corresponding to `x^r` and the determinent of the matrix A =
350 [A_ij]. In other words, all of the fixed (cachable) data needed
351 to compute the coefficients of the characteristic polynomial.
356 # Construct a new algebra over a multivariate polynomial ring...
357 names
= ['X' + str(i
) for i
in range(1,n
+1)]
358 R
= PolynomialRing(self
.base_ring(), names
)
359 J
= FiniteDimensionalEuclideanJordanAlgebra(R
,
360 self
._multiplication
_table
,
363 idmat
= identity_matrix(J
.base_ring(), n
)
365 W
= self
._charpoly
_basis
_space
()
366 W
= W
.change_ring(R
.fraction_field())
368 # Starting with the standard coordinates x = (X1,X2,...,Xn)
369 # and then converting the entries to W-coordinates allows us
370 # to pass in the standard coordinates to the charpoly and get
371 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
374 # W.coordinates(x^2) eval'd at (standard z-coords)
378 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
380 # We want the middle equivalent thing in our matrix, but use
381 # the first equivalent thing instead so that we can pass in
382 # standard coordinates.
383 x
= J(vector(R
, R
.gens()))
384 l1
= [column_matrix(W
.coordinates((x
**k
).vector())) for k
in range(r
)]
385 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
386 A_of_x
= block_matrix(R
, 1, n
, (l1
+ l2
))
387 xr
= W
.coordinates((x
**r
).vector())
388 return (A_of_x
, x
, xr
, A_of_x
.det())
392 def characteristic_polynomial(self
):
397 This implementation doesn't guarantee that the polynomial
398 denominator in the coefficients is not identically zero, so
399 theoretically it could crash. The way that this is handled
400 in e.g. Faraut and Koranyi is to use a basis that guarantees
401 the denominator is non-zero. But, doing so requires knowledge
402 of at least one regular element, and we don't even know how
403 to do that. The trade-off is that, if we use the standard basis,
404 the resulting polynomial will accept the "usual" coordinates. In
405 other words, we don't have to do a change of basis before e.g.
406 computing the trace or determinant.
410 The characteristic polynomial in the spin algebra is given in
411 Alizadeh, Example 11.11::
413 sage: J = JordanSpinEJA(3)
414 sage: p = J.characteristic_polynomial(); p
415 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
416 sage: xvec = J.one().vector()
424 # The list of coefficient polynomials a_1, a_2, ..., a_n.
425 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
427 # We go to a bit of trouble here to reorder the
428 # indeterminates, so that it's easier to evaluate the
429 # characteristic polynomial at x's coordinates and get back
430 # something in terms of t, which is what we want.
432 S
= PolynomialRing(self
.base_ring(),'t')
434 S
= PolynomialRing(S
, R
.variable_names())
437 # Note: all entries past the rth should be zero. The
438 # coefficient of the highest power (x^r) is 1, but it doesn't
439 # appear in the solution vector which contains coefficients
440 # for the other powers (to make them sum to x^r).
442 a
[r
] = 1 # corresponds to x^r
444 # When the rank is equal to the dimension, trying to
445 # assign a[r] goes out-of-bounds.
446 a
.append(1) # corresponds to x^r
448 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
451 def inner_product(self
, x
, y
):
453 The inner product associated with this Euclidean Jordan algebra.
455 Defaults to the trace inner product, but can be overridden by
456 subclasses if they are sure that the necessary properties are
461 The inner product must satisfy its axiom for this algebra to truly
462 be a Euclidean Jordan Algebra::
464 sage: set_random_seed()
465 sage: J = random_eja()
466 sage: x = J.random_element()
467 sage: y = J.random_element()
468 sage: z = J.random_element()
469 sage: (x*y).inner_product(z) == y.inner_product(x*z)
473 if (not x
in self
) or (not y
in self
):
474 raise TypeError("arguments must live in this algebra")
475 return x
.trace_inner_product(y
)
478 def natural_basis(self
):
480 Return a more-natural representation of this algebra's basis.
482 Every finite-dimensional Euclidean Jordan Algebra is a direct
483 sum of five simple algebras, four of which comprise Hermitian
484 matrices. This method returns the original "natural" basis
485 for our underlying vector space. (Typically, the natural basis
486 is used to construct the multiplication table in the first place.)
488 Note that this will always return a matrix. The standard basis
489 in `R^n` will be returned as `n`-by-`1` column matrices.
493 sage: J = RealSymmetricEJA(2)
496 sage: J.natural_basis()
504 sage: J = JordanSpinEJA(2)
507 sage: J.natural_basis()
514 if self
._natural
_basis
is None:
515 return tuple( b
.vector().column() for b
in self
.basis() )
517 return self
._natural
_basis
522 Return the rank of this EJA.
524 if self
._rank
is None:
525 raise ValueError("no rank specified at genesis")
530 class Element(FiniteDimensionalAlgebraElement
):
532 An element of a Euclidean Jordan algebra.
537 Oh man, I should not be doing this. This hides the "disabled"
538 methods ``left_matrix`` and ``matrix`` from introspection;
539 in particular it removes them from tab-completion.
541 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
542 dir(self
.__class
__) )
545 def __init__(self
, A
, elt
=None):
549 The identity in `S^n` is converted to the identity in the EJA::
551 sage: J = RealSymmetricEJA(3)
552 sage: I = identity_matrix(QQ,3)
553 sage: J(I) == J.one()
556 This skew-symmetric matrix can't be represented in the EJA::
558 sage: J = RealSymmetricEJA(3)
559 sage: A = matrix(QQ,3, lambda i,j: i-j)
561 Traceback (most recent call last):
563 ArithmeticError: vector is not in free module
566 # Goal: if we're given a matrix, and if it lives in our
567 # parent algebra's "natural ambient space," convert it
568 # into an algebra element.
570 # The catch is, we make a recursive call after converting
571 # the given matrix into a vector that lives in the algebra.
572 # This we need to try the parent class initializer first,
573 # to avoid recursing forever if we're given something that
574 # already fits into the algebra, but also happens to live
575 # in the parent's "natural ambient space" (this happens with
578 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, elt
)
580 natural_basis
= A
.natural_basis()
581 if elt
in natural_basis
[0].matrix_space():
582 # Thanks for nothing! Matrix spaces aren't vector
583 # spaces in Sage, so we have to figure out its
584 # natural-basis coordinates ourselves.
585 V
= VectorSpace(elt
.base_ring(), elt
.nrows()**2)
586 W
= V
.span( _mat2vec(s
) for s
in natural_basis
)
587 coords
= W
.coordinates(_mat2vec(elt
))
588 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, coords
)
590 def __pow__(self
, n
):
592 Return ``self`` raised to the power ``n``.
594 Jordan algebras are always power-associative; see for
595 example Faraut and Koranyi, Proposition II.1.2 (ii).
599 We have to override this because our superclass uses row vectors
600 instead of column vectors! We, on the other hand, assume column
605 sage: set_random_seed()
606 sage: x = random_eja().random_element()
607 sage: x.operator_matrix()*x.vector() == (x^2).vector()
610 A few examples of power-associativity::
612 sage: set_random_seed()
613 sage: x = random_eja().random_element()
614 sage: x*(x*x)*(x*x) == x^5
616 sage: (x*x)*(x*x*x) == x^5
619 We also know that powers operator-commute (Koecher, Chapter
622 sage: set_random_seed()
623 sage: x = random_eja().random_element()
624 sage: m = ZZ.random_element(0,10)
625 sage: n = ZZ.random_element(0,10)
626 sage: Lxm = (x^m).operator_matrix()
627 sage: Lxn = (x^n).operator_matrix()
628 sage: Lxm*Lxn == Lxn*Lxm
638 return A( (self
.operator_matrix()**(n
-1))*self
.vector() )
641 def apply_univariate_polynomial(self
, p
):
643 Apply the univariate polynomial ``p`` to this element.
645 A priori, SageMath won't allow us to apply a univariate
646 polynomial to an element of an EJA, because we don't know
647 that EJAs are rings (they are usually not associative). Of
648 course, we know that EJAs are power-associative, so the
649 operation is ultimately kosher. This function sidesteps
650 the CAS to get the answer we want and expect.
654 sage: R = PolynomialRing(QQ, 't')
656 sage: p = t^4 - t^3 + 5*t - 2
657 sage: J = RealCartesianProductEJA(5)
658 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
663 We should always get back an element of the algebra::
665 sage: set_random_seed()
666 sage: p = PolynomialRing(QQ, 't').random_element()
667 sage: J = random_eja()
668 sage: x = J.random_element()
669 sage: x.apply_univariate_polynomial(p) in J
673 if len(p
.variables()) > 1:
674 raise ValueError("not a univariate polynomial")
677 # Convert the coeficcients to the parent's base ring,
678 # because a priori they might live in an (unnecessarily)
679 # larger ring for which P.sum() would fail below.
680 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
681 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
684 def characteristic_polynomial(self
):
686 Return the characteristic polynomial of this element.
690 The rank of `R^3` is three, and the minimal polynomial of
691 the identity element is `(t-1)` from which it follows that
692 the characteristic polynomial should be `(t-1)^3`::
694 sage: J = RealCartesianProductEJA(3)
695 sage: J.one().characteristic_polynomial()
696 t^3 - 3*t^2 + 3*t - 1
698 Likewise, the characteristic of the zero element in the
699 rank-three algebra `R^{n}` should be `t^{3}`::
701 sage: J = RealCartesianProductEJA(3)
702 sage: J.zero().characteristic_polynomial()
705 The characteristic polynomial of an element should evaluate
706 to zero on that element::
708 sage: set_random_seed()
709 sage: x = RealCartesianProductEJA(3).random_element()
710 sage: p = x.characteristic_polynomial()
711 sage: x.apply_univariate_polynomial(p)
715 p
= self
.parent().characteristic_polynomial()
716 return p(*self
.vector())
719 def inner_product(self
, other
):
721 Return the parent algebra's inner product of myself and ``other``.
725 The inner product in the Jordan spin algebra is the usual
726 inner product on `R^n` (this example only works because the
727 basis for the Jordan algebra is the standard basis in `R^n`)::
729 sage: J = JordanSpinEJA(3)
730 sage: x = vector(QQ,[1,2,3])
731 sage: y = vector(QQ,[4,5,6])
732 sage: x.inner_product(y)
734 sage: J(x).inner_product(J(y))
737 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
738 multiplication is the usual matrix multiplication in `S^n`,
739 so the inner product of the identity matrix with itself
742 sage: J = RealSymmetricEJA(3)
743 sage: J.one().inner_product(J.one())
746 Likewise, the inner product on `C^n` is `<X,Y> =
747 Re(trace(X*Y))`, where we must necessarily take the real
748 part because the product of Hermitian matrices may not be
751 sage: J = ComplexHermitianEJA(3)
752 sage: J.one().inner_product(J.one())
755 Ditto for the quaternions::
757 sage: J = QuaternionHermitianEJA(3)
758 sage: J.one().inner_product(J.one())
763 Ensure that we can always compute an inner product, and that
764 it gives us back a real number::
766 sage: set_random_seed()
767 sage: J = random_eja()
768 sage: x = J.random_element()
769 sage: y = J.random_element()
770 sage: x.inner_product(y) in RR
776 raise TypeError("'other' must live in the same algebra")
778 return P
.inner_product(self
, other
)
781 def operator_commutes_with(self
, other
):
783 Return whether or not this element operator-commutes
788 The definition of a Jordan algebra says that any element
789 operator-commutes with its square::
791 sage: set_random_seed()
792 sage: x = random_eja().random_element()
793 sage: x.operator_commutes_with(x^2)
798 Test Lemma 1 from Chapter III of Koecher::
800 sage: set_random_seed()
801 sage: J = random_eja()
802 sage: u = J.random_element()
803 sage: v = J.random_element()
804 sage: lhs = u.operator_commutes_with(u*v)
805 sage: rhs = v.operator_commutes_with(u^2)
810 if not other
in self
.parent():
811 raise TypeError("'other' must live in the same algebra")
813 A
= self
.operator_matrix()
814 B
= other
.operator_matrix()
820 Return my determinant, the product of my eigenvalues.
824 sage: J = JordanSpinEJA(2)
825 sage: e0,e1 = J.gens()
826 sage: x = sum( J.gens() )
832 sage: J = JordanSpinEJA(3)
833 sage: e0,e1,e2 = J.gens()
834 sage: x = sum( J.gens() )
840 An element is invertible if and only if its determinant is
843 sage: set_random_seed()
844 sage: x = random_eja().random_element()
845 sage: x.is_invertible() == (x.det() != 0)
851 p
= P
._charpoly
_coeff
(0)
852 # The _charpoly_coeff function already adds the factor of
853 # -1 to ensure that _charpoly_coeff(0) is really what
854 # appears in front of t^{0} in the charpoly. However,
855 # we want (-1)^r times THAT for the determinant.
856 return ((-1)**r
)*p(*self
.vector())
861 Return the Jordan-multiplicative inverse of this element.
865 We appeal to the quadratic representation as in Koecher's
866 Theorem 12 in Chapter III, Section 5.
870 The inverse in the spin factor algebra is given in Alizadeh's
873 sage: set_random_seed()
874 sage: n = ZZ.random_element(1,10)
875 sage: J = JordanSpinEJA(n)
876 sage: x = J.random_element()
877 sage: while x.is_zero():
878 ....: x = J.random_element()
879 sage: x_vec = x.vector()
881 sage: x_bar = x_vec[1:]
882 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
883 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
884 sage: x_inverse = coeff*inv_vec
885 sage: x.inverse() == J(x_inverse)
890 The identity element is its own inverse::
892 sage: set_random_seed()
893 sage: J = random_eja()
894 sage: J.one().inverse() == J.one()
897 If an element has an inverse, it acts like one::
899 sage: set_random_seed()
900 sage: J = random_eja()
901 sage: x = J.random_element()
902 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
905 The inverse of the inverse is what we started with::
907 sage: set_random_seed()
908 sage: J = random_eja()
909 sage: x = J.random_element()
910 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
913 The zero element is never invertible::
915 sage: set_random_seed()
916 sage: J = random_eja().zero().inverse()
917 Traceback (most recent call last):
919 ValueError: element is not invertible
922 if not self
.is_invertible():
923 raise ValueError("element is not invertible")
926 return P(self
.quadratic_representation().inverse()*self
.vector())
929 def is_invertible(self
):
931 Return whether or not this element is invertible.
933 We can't use the superclass method because it relies on
934 the algebra being associative.
938 The usual way to do this is to check if the determinant is
939 zero, but we need the characteristic polynomial for the
940 determinant. The minimal polynomial is a lot easier to get,
941 so we use Corollary 2 in Chapter V of Koecher to check
942 whether or not the paren't algebra's zero element is a root
943 of this element's minimal polynomial.
947 The identity element is always invertible::
949 sage: set_random_seed()
950 sage: J = random_eja()
951 sage: J.one().is_invertible()
954 The zero element is never invertible::
956 sage: set_random_seed()
957 sage: J = random_eja()
958 sage: J.zero().is_invertible()
962 zero
= self
.parent().zero()
963 p
= self
.minimal_polynomial()
964 return not (p(zero
) == zero
)
967 def is_nilpotent(self
):
969 Return whether or not some power of this element is zero.
971 The superclass method won't work unless we're in an
972 associative algebra, and we aren't. However, we generate
973 an assocoative subalgebra and we're nilpotent there if and
974 only if we're nilpotent here (probably).
978 The identity element is never nilpotent::
980 sage: set_random_seed()
981 sage: random_eja().one().is_nilpotent()
984 The additive identity is always nilpotent::
986 sage: set_random_seed()
987 sage: random_eja().zero().is_nilpotent()
991 # The element we're going to call "is_nilpotent()" on.
992 # Either myself, interpreted as an element of a finite-
993 # dimensional algebra, or an element of an associative
997 if self
.parent().is_associative():
998 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1000 V
= self
.span_of_powers()
1001 assoc_subalg
= self
.subalgebra_generated_by()
1002 # Mis-design warning: the basis used for span_of_powers()
1003 # and subalgebra_generated_by() must be the same, and in
1005 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1007 # Recursive call, but should work since elt lives in an
1008 # associative algebra.
1009 return elt
.is_nilpotent()
1012 def is_regular(self
):
1014 Return whether or not this is a regular element.
1018 The identity element always has degree one, but any element
1019 linearly-independent from it is regular::
1021 sage: J = JordanSpinEJA(5)
1022 sage: J.one().is_regular()
1024 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
1025 sage: for x in J.gens():
1026 ....: (J.one() + x).is_regular()
1034 return self
.degree() == self
.parent().rank()
1039 Compute the degree of this element the straightforward way
1040 according to the definition; by appending powers to a list
1041 and figuring out its dimension (that is, whether or not
1042 they're linearly dependent).
1046 sage: J = JordanSpinEJA(4)
1047 sage: J.one().degree()
1049 sage: e0,e1,e2,e3 = J.gens()
1050 sage: (e0 - e1).degree()
1053 In the spin factor algebra (of rank two), all elements that
1054 aren't multiples of the identity are regular::
1056 sage: set_random_seed()
1057 sage: n = ZZ.random_element(1,10)
1058 sage: J = JordanSpinEJA(n)
1059 sage: x = J.random_element()
1060 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
1064 return self
.span_of_powers().dimension()
1067 def left_matrix(self
):
1069 Our parent class defines ``left_matrix`` and ``matrix``
1070 methods whose names are misleading. We don't want them.
1072 raise NotImplementedError("use operator_matrix() instead")
1074 matrix
= left_matrix
1077 def minimal_polynomial(self
):
1079 Return the minimal polynomial of this element,
1080 as a function of the variable `t`.
1084 We restrict ourselves to the associative subalgebra
1085 generated by this element, and then return the minimal
1086 polynomial of this element's operator matrix (in that
1087 subalgebra). This works by Baes Proposition 2.3.16.
1091 The minimal polynomial of the identity and zero elements are
1094 sage: set_random_seed()
1095 sage: J = random_eja()
1096 sage: J.one().minimal_polynomial()
1098 sage: J.zero().minimal_polynomial()
1101 The degree of an element is (by one definition) the degree
1102 of its minimal polynomial::
1104 sage: set_random_seed()
1105 sage: x = random_eja().random_element()
1106 sage: x.degree() == x.minimal_polynomial().degree()
1109 The minimal polynomial and the characteristic polynomial coincide
1110 and are known (see Alizadeh, Example 11.11) for all elements of
1111 the spin factor algebra that aren't scalar multiples of the
1114 sage: set_random_seed()
1115 sage: n = ZZ.random_element(2,10)
1116 sage: J = JordanSpinEJA(n)
1117 sage: y = J.random_element()
1118 sage: while y == y.coefficient(0)*J.one():
1119 ....: y = J.random_element()
1120 sage: y0 = y.vector()[0]
1121 sage: y_bar = y.vector()[1:]
1122 sage: actual = y.minimal_polynomial()
1123 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1124 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1125 sage: bool(actual == expected)
1128 The minimal polynomial should always kill its element::
1130 sage: set_random_seed()
1131 sage: x = random_eja().random_element()
1132 sage: p = x.minimal_polynomial()
1133 sage: x.apply_univariate_polynomial(p)
1137 V
= self
.span_of_powers()
1138 assoc_subalg
= self
.subalgebra_generated_by()
1139 # Mis-design warning: the basis used for span_of_powers()
1140 # and subalgebra_generated_by() must be the same, and in
1142 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1144 # We get back a symbolic polynomial in 'x' but want a real
1145 # polynomial in 't'.
1146 p_of_x
= elt
.operator_matrix().minimal_polynomial()
1147 return p_of_x
.change_variable_name('t')
1150 def natural_representation(self
):
1152 Return a more-natural representation of this element.
1154 Every finite-dimensional Euclidean Jordan Algebra is a
1155 direct sum of five simple algebras, four of which comprise
1156 Hermitian matrices. This method returns the original
1157 "natural" representation of this element as a Hermitian
1158 matrix, if it has one. If not, you get the usual representation.
1162 sage: J = ComplexHermitianEJA(3)
1165 sage: J.one().natural_representation()
1175 sage: J = QuaternionHermitianEJA(3)
1178 sage: J.one().natural_representation()
1179 [1 0 0 0 0 0 0 0 0 0 0 0]
1180 [0 1 0 0 0 0 0 0 0 0 0 0]
1181 [0 0 1 0 0 0 0 0 0 0 0 0]
1182 [0 0 0 1 0 0 0 0 0 0 0 0]
1183 [0 0 0 0 1 0 0 0 0 0 0 0]
1184 [0 0 0 0 0 1 0 0 0 0 0 0]
1185 [0 0 0 0 0 0 1 0 0 0 0 0]
1186 [0 0 0 0 0 0 0 1 0 0 0 0]
1187 [0 0 0 0 0 0 0 0 1 0 0 0]
1188 [0 0 0 0 0 0 0 0 0 1 0 0]
1189 [0 0 0 0 0 0 0 0 0 0 1 0]
1190 [0 0 0 0 0 0 0 0 0 0 0 1]
1193 B
= self
.parent().natural_basis()
1194 W
= B
[0].matrix_space()
1195 return W
.linear_combination(zip(self
.vector(), B
))
1200 Return the left-multiplication-by-this-element
1201 operator on the ambient algebra.
1205 sage: set_random_seed()
1206 sage: J = random_eja()
1207 sage: x = J.random_element()
1208 sage: y = J.random_element()
1209 sage: x.operator()(y) == x*y
1211 sage: y.operator()(x) == x*y
1216 return FiniteDimensionalEuclideanJordanAlgebraMorphism(
1218 self
.operator_matrix() )
1222 def operator_matrix(self
):
1224 Return the matrix that represents left- (or right-)
1225 multiplication by this element in the parent algebra.
1227 We implement this ourselves to work around the fact that
1228 our parent class represents everything with row vectors.
1232 Test the first polarization identity from my notes, Koecher Chapter
1233 III, or from Baes (2.3)::
1235 sage: set_random_seed()
1236 sage: J = random_eja()
1237 sage: x = J.random_element()
1238 sage: y = J.random_element()
1239 sage: Lx = x.operator_matrix()
1240 sage: Ly = y.operator_matrix()
1241 sage: Lxx = (x*x).operator_matrix()
1242 sage: Lxy = (x*y).operator_matrix()
1243 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
1246 Test the second polarization identity from my notes or from
1249 sage: set_random_seed()
1250 sage: J = random_eja()
1251 sage: x = J.random_element()
1252 sage: y = J.random_element()
1253 sage: z = J.random_element()
1254 sage: Lx = x.operator_matrix()
1255 sage: Ly = y.operator_matrix()
1256 sage: Lz = z.operator_matrix()
1257 sage: Lzy = (z*y).operator_matrix()
1258 sage: Lxy = (x*y).operator_matrix()
1259 sage: Lxz = (x*z).operator_matrix()
1260 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
1263 Test the third polarization identity from my notes or from
1266 sage: set_random_seed()
1267 sage: J = random_eja()
1268 sage: u = J.random_element()
1269 sage: y = J.random_element()
1270 sage: z = J.random_element()
1271 sage: Lu = u.operator_matrix()
1272 sage: Ly = y.operator_matrix()
1273 sage: Lz = z.operator_matrix()
1274 sage: Lzy = (z*y).operator_matrix()
1275 sage: Luy = (u*y).operator_matrix()
1276 sage: Luz = (u*z).operator_matrix()
1277 sage: Luyz = (u*(y*z)).operator_matrix()
1278 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
1279 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
1280 sage: bool(lhs == rhs)
1284 fda_elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1285 return fda_elt
.matrix().transpose()
1288 def quadratic_representation(self
, other
=None):
1290 Return the quadratic representation of this element.
1294 The explicit form in the spin factor algebra is given by
1295 Alizadeh's Example 11.12::
1297 sage: set_random_seed()
1298 sage: n = ZZ.random_element(1,10)
1299 sage: J = JordanSpinEJA(n)
1300 sage: x = J.random_element()
1301 sage: x_vec = x.vector()
1303 sage: x_bar = x_vec[1:]
1304 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1305 sage: B = 2*x0*x_bar.row()
1306 sage: C = 2*x0*x_bar.column()
1307 sage: D = identity_matrix(QQ, n-1)
1308 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1309 sage: D = D + 2*x_bar.tensor_product(x_bar)
1310 sage: Q = block_matrix(2,2,[A,B,C,D])
1311 sage: Q == x.quadratic_representation()
1314 Test all of the properties from Theorem 11.2 in Alizadeh::
1316 sage: set_random_seed()
1317 sage: J = random_eja()
1318 sage: x = J.random_element()
1319 sage: y = J.random_element()
1320 sage: Lx = x.operator_matrix()
1321 sage: Lxx = (x*x).operator_matrix()
1322 sage: Qx = x.quadratic_representation()
1323 sage: Qy = y.quadratic_representation()
1324 sage: Qxy = x.quadratic_representation(y)
1325 sage: Qex = J.one().quadratic_representation(x)
1326 sage: n = ZZ.random_element(10)
1327 sage: Qxn = (x^n).quadratic_representation()
1331 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1336 sage: alpha = QQ.random_element()
1337 sage: (alpha*x).quadratic_representation() == (alpha^2)*Qx
1342 sage: not x.is_invertible() or (
1343 ....: Qx*x.inverse().vector() == x.vector() )
1346 sage: not x.is_invertible() or (
1349 ....: x.inverse().quadratic_representation() )
1352 sage: Qxy*(J.one().vector()) == (x*y).vector()
1357 sage: not x.is_invertible() or (
1358 ....: x.quadratic_representation(x.inverse())*Qx
1359 ....: == Qx*x.quadratic_representation(x.inverse()) )
1362 sage: not x.is_invertible() or (
1363 ....: x.quadratic_representation(x.inverse())*Qx
1365 ....: 2*x.operator_matrix()*Qex - Qx )
1368 sage: 2*x.operator_matrix()*Qex - Qx == Lxx
1373 sage: J(Qy*x.vector()).quadratic_representation() == Qy*Qx*Qy
1383 sage: not x.is_invertible() or (
1384 ....: Qx*x.inverse().operator_matrix() == Lx )
1389 sage: not x.operator_commutes_with(y) or (
1390 ....: J(Qx*y.vector())^n == J(Qxn*(y^n).vector()) )
1396 elif not other
in self
.parent():
1397 raise TypeError("'other' must live in the same algebra")
1399 L
= self
.operator_matrix()
1400 M
= other
.operator_matrix()
1401 return ( L
*M
+ M
*L
- (self
*other
).operator_matrix() )
1404 def span_of_powers(self
):
1406 Return the vector space spanned by successive powers of
1409 # The dimension of the subalgebra can't be greater than
1410 # the big algebra, so just put everything into a list
1411 # and let span() get rid of the excess.
1413 # We do the extra ambient_vector_space() in case we're messing
1414 # with polynomials and the direct parent is a module.
1415 V
= self
.vector().parent().ambient_vector_space()
1416 return V
.span( (self
**d
).vector() for d
in xrange(V
.dimension()) )
1419 def subalgebra_generated_by(self
):
1421 Return the associative subalgebra of the parent EJA generated
1426 sage: set_random_seed()
1427 sage: x = random_eja().random_element()
1428 sage: x.subalgebra_generated_by().is_associative()
1431 Squaring in the subalgebra should be the same thing as
1432 squaring in the superalgebra::
1434 sage: set_random_seed()
1435 sage: x = random_eja().random_element()
1436 sage: u = x.subalgebra_generated_by().random_element()
1437 sage: u.operator_matrix()*u.vector() == (u**2).vector()
1441 # First get the subspace spanned by the powers of myself...
1442 V
= self
.span_of_powers()
1443 F
= self
.base_ring()
1445 # Now figure out the entries of the right-multiplication
1446 # matrix for the successive basis elements b0, b1,... of
1449 for b_right
in V
.basis():
1450 eja_b_right
= self
.parent()(b_right
)
1452 # The first row of the right-multiplication matrix by
1453 # b1 is what we get if we apply that matrix to b1. The
1454 # second row of the right multiplication matrix by b1
1455 # is what we get when we apply that matrix to b2...
1457 # IMPORTANT: this assumes that all vectors are COLUMN
1458 # vectors, unlike our superclass (which uses row vectors).
1459 for b_left
in V
.basis():
1460 eja_b_left
= self
.parent()(b_left
)
1461 # Multiply in the original EJA, but then get the
1462 # coordinates from the subalgebra in terms of its
1464 this_row
= V
.coordinates((eja_b_left
*eja_b_right
).vector())
1465 b_right_rows
.append(this_row
)
1466 b_right_matrix
= matrix(F
, b_right_rows
)
1467 mats
.append(b_right_matrix
)
1469 # It's an algebra of polynomials in one element, and EJAs
1470 # are power-associative.
1472 # TODO: choose generator names intelligently.
1473 return FiniteDimensionalEuclideanJordanAlgebra(F
, mats
, assume_associative
=True, names
='f')
1476 def subalgebra_idempotent(self
):
1478 Find an idempotent in the associative subalgebra I generate
1479 using Proposition 2.3.5 in Baes.
1483 sage: set_random_seed()
1484 sage: J = random_eja()
1485 sage: x = J.random_element()
1486 sage: while x.is_nilpotent():
1487 ....: x = J.random_element()
1488 sage: c = x.subalgebra_idempotent()
1493 if self
.is_nilpotent():
1494 raise ValueError("this only works with non-nilpotent elements!")
1496 V
= self
.span_of_powers()
1497 J
= self
.subalgebra_generated_by()
1498 # Mis-design warning: the basis used for span_of_powers()
1499 # and subalgebra_generated_by() must be the same, and in
1501 u
= J(V
.coordinates(self
.vector()))
1503 # The image of the matrix of left-u^m-multiplication
1504 # will be minimal for some natural number s...
1506 minimal_dim
= V
.dimension()
1507 for i
in xrange(1, V
.dimension()):
1508 this_dim
= (u
**i
).operator_matrix().image().dimension()
1509 if this_dim
< minimal_dim
:
1510 minimal_dim
= this_dim
1513 # Now minimal_matrix should correspond to the smallest
1514 # non-zero subspace in Baes's (or really, Koecher's)
1517 # However, we need to restrict the matrix to work on the
1518 # subspace... or do we? Can't we just solve, knowing that
1519 # A(c) = u^(s+1) should have a solution in the big space,
1522 # Beware, solve_right() means that we're using COLUMN vectors.
1523 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1525 A
= u_next
.operator_matrix()
1526 c_coordinates
= A
.solve_right(u_next
.vector())
1528 # Now c_coordinates is the idempotent we want, but it's in
1529 # the coordinate system of the subalgebra.
1531 # We need the basis for J, but as elements of the parent algebra.
1533 basis
= [self
.parent(v
) for v
in V
.basis()]
1534 return self
.parent().linear_combination(zip(c_coordinates
, basis
))
1539 Return my trace, the sum of my eigenvalues.
1543 sage: J = JordanSpinEJA(3)
1544 sage: x = sum(J.gens())
1550 sage: J = RealCartesianProductEJA(5)
1551 sage: J.one().trace()
1556 The trace of an element is a real number::
1558 sage: set_random_seed()
1559 sage: J = random_eja()
1560 sage: J.random_element().trace() in J.base_ring()
1566 p
= P
._charpoly
_coeff
(r
-1)
1567 # The _charpoly_coeff function already adds the factor of
1568 # -1 to ensure that _charpoly_coeff(r-1) is really what
1569 # appears in front of t^{r-1} in the charpoly. However,
1570 # we want the negative of THAT for the trace.
1571 return -p(*self
.vector())
1574 def trace_inner_product(self
, other
):
1576 Return the trace inner product of myself and ``other``.
1580 The trace inner product is commutative::
1582 sage: set_random_seed()
1583 sage: J = random_eja()
1584 sage: x = J.random_element(); y = J.random_element()
1585 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1588 The trace inner product is bilinear::
1590 sage: set_random_seed()
1591 sage: J = random_eja()
1592 sage: x = J.random_element()
1593 sage: y = J.random_element()
1594 sage: z = J.random_element()
1595 sage: a = QQ.random_element();
1596 sage: actual = (a*(x+z)).trace_inner_product(y)
1597 sage: expected = ( a*x.trace_inner_product(y) +
1598 ....: a*z.trace_inner_product(y) )
1599 sage: actual == expected
1601 sage: actual = x.trace_inner_product(a*(y+z))
1602 sage: expected = ( a*x.trace_inner_product(y) +
1603 ....: a*x.trace_inner_product(z) )
1604 sage: actual == expected
1607 The trace inner product satisfies the compatibility
1608 condition in the definition of a Euclidean Jordan algebra::
1610 sage: set_random_seed()
1611 sage: J = random_eja()
1612 sage: x = J.random_element()
1613 sage: y = J.random_element()
1614 sage: z = J.random_element()
1615 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1619 if not other
in self
.parent():
1620 raise TypeError("'other' must live in the same algebra")
1622 return (self
*other
).trace()
1625 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1627 Return the Euclidean Jordan Algebra corresponding to the set
1628 `R^n` under the Hadamard product.
1630 Note: this is nothing more than the Cartesian product of ``n``
1631 copies of the spin algebra. Once Cartesian product algebras
1632 are implemented, this can go.
1636 This multiplication table can be verified by hand::
1638 sage: J = RealCartesianProductEJA(3)
1639 sage: e0,e1,e2 = J.gens()
1655 def __classcall_private__(cls
, n
, field
=QQ
):
1656 # The FiniteDimensionalAlgebra constructor takes a list of
1657 # matrices, the ith representing right multiplication by the ith
1658 # basis element in the vector space. So if e_1 = (1,0,0), then
1659 # right (Hadamard) multiplication of x by e_1 picks out the first
1660 # component of x; and likewise for the ith basis element e_i.
1661 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
1662 for i
in xrange(n
) ]
1664 fdeja
= super(RealCartesianProductEJA
, cls
)
1665 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=n
)
1667 def inner_product(self
, x
, y
):
1668 return _usual_ip(x
,y
)
1673 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1677 For now, we choose a random natural number ``n`` (greater than zero)
1678 and then give you back one of the following:
1680 * The cartesian product of the rational numbers ``n`` times; this is
1681 ``QQ^n`` with the Hadamard product.
1683 * The Jordan spin algebra on ``QQ^n``.
1685 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1688 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1689 in the space of ``2n``-by-``2n`` real symmetric matrices.
1691 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1692 in the space of ``4n``-by-``4n`` real symmetric matrices.
1694 Later this might be extended to return Cartesian products of the
1700 Euclidean Jordan algebra of degree...
1704 # The max_n component lets us choose different upper bounds on the
1705 # value "n" that gets passed to the constructor. This is needed
1706 # because e.g. R^{10} is reasonable to test, while the Hermitian
1707 # 10-by-10 quaternion matrices are not.
1708 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
1710 (RealSymmetricEJA
, 5),
1711 (ComplexHermitianEJA
, 4),
1712 (QuaternionHermitianEJA
, 3)])
1713 n
= ZZ
.random_element(1, max_n
)
1714 return constructor(n
, field
=QQ
)
1718 def _real_symmetric_basis(n
, field
=QQ
):
1720 Return a basis for the space of real symmetric n-by-n matrices.
1722 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1726 for j
in xrange(i
+1):
1727 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1731 # Beware, orthogonal but not normalized!
1732 Sij
= Eij
+ Eij
.transpose()
1737 def _complex_hermitian_basis(n
, field
=QQ
):
1739 Returns a basis for the space of complex Hermitian n-by-n matrices.
1743 sage: set_random_seed()
1744 sage: n = ZZ.random_element(1,5)
1745 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1749 F
= QuadraticField(-1, 'I')
1752 # This is like the symmetric case, but we need to be careful:
1754 # * We want conjugate-symmetry, not just symmetry.
1755 # * The diagonal will (as a result) be real.
1759 for j
in xrange(i
+1):
1760 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1762 Sij
= _embed_complex_matrix(Eij
)
1765 # Beware, orthogonal but not normalized! The second one
1766 # has a minus because it's conjugated.
1767 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
1769 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
1774 def _quaternion_hermitian_basis(n
, field
=QQ
):
1776 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1780 sage: set_random_seed()
1781 sage: n = ZZ.random_element(1,5)
1782 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1786 Q
= QuaternionAlgebra(QQ
,-1,-1)
1789 # This is like the symmetric case, but we need to be careful:
1791 # * We want conjugate-symmetry, not just symmetry.
1792 # * The diagonal will (as a result) be real.
1796 for j
in xrange(i
+1):
1797 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1799 Sij
= _embed_quaternion_matrix(Eij
)
1802 # Beware, orthogonal but not normalized! The second,
1803 # third, and fourth ones have a minus because they're
1805 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
1807 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
1809 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
1811 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
1817 return vector(m
.base_ring(), m
.list())
1820 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
1822 def _multiplication_table_from_matrix_basis(basis
):
1824 At least three of the five simple Euclidean Jordan algebras have the
1825 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1826 multiplication on the right is matrix multiplication. Given a basis
1827 for the underlying matrix space, this function returns a
1828 multiplication table (obtained by looping through the basis
1829 elements) for an algebra of those matrices. A reordered copy
1830 of the basis is also returned to work around the fact that
1831 the ``span()`` in this function will change the order of the basis
1832 from what we think it is, to... something else.
1834 # In S^2, for example, we nominally have four coordinates even
1835 # though the space is of dimension three only. The vector space V
1836 # is supposed to hold the entire long vector, and the subspace W
1837 # of V will be spanned by the vectors that arise from symmetric
1838 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1839 field
= basis
[0].base_ring()
1840 dimension
= basis
[0].nrows()
1842 V
= VectorSpace(field
, dimension
**2)
1843 W
= V
.span( _mat2vec(s
) for s
in basis
)
1845 # Taking the span above reorders our basis (thanks, jerk!) so we
1846 # need to put our "matrix basis" in the same order as the
1847 # (reordered) vector basis.
1848 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
1852 # Brute force the multiplication-by-s matrix by looping
1853 # through all elements of the basis and doing the computation
1854 # to find out what the corresponding row should be. BEWARE:
1855 # these multiplication tables won't be symmetric! It therefore
1856 # becomes REALLY IMPORTANT that the underlying algebra
1857 # constructor uses ROW vectors and not COLUMN vectors. That's
1858 # why we're computing rows here and not columns.
1861 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
1862 Q_rows
.append(W
.coordinates(this_row
))
1863 Q
= matrix(field
, W
.dimension(), Q_rows
)
1869 def _embed_complex_matrix(M
):
1871 Embed the n-by-n complex matrix ``M`` into the space of real
1872 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1873 bi` to the block matrix ``[[a,b],[-b,a]]``.
1877 sage: F = QuadraticField(-1,'i')
1878 sage: x1 = F(4 - 2*i)
1879 sage: x2 = F(1 + 2*i)
1882 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1883 sage: _embed_complex_matrix(M)
1892 Embedding is a homomorphism (isomorphism, in fact)::
1894 sage: set_random_seed()
1895 sage: n = ZZ.random_element(5)
1896 sage: F = QuadraticField(-1, 'i')
1897 sage: X = random_matrix(F, n)
1898 sage: Y = random_matrix(F, n)
1899 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1900 sage: expected = _embed_complex_matrix(X*Y)
1901 sage: actual == expected
1907 raise ValueError("the matrix 'M' must be square")
1908 field
= M
.base_ring()
1913 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1915 # We can drop the imaginaries here.
1916 return block_matrix(field
.base_ring(), n
, blocks
)
1919 def _unembed_complex_matrix(M
):
1921 The inverse of _embed_complex_matrix().
1925 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1926 ....: [-2, 1, -4, 3],
1927 ....: [ 9, 10, 11, 12],
1928 ....: [-10, 9, -12, 11] ])
1929 sage: _unembed_complex_matrix(A)
1931 [ 10*i + 9 12*i + 11]
1935 Unembedding is the inverse of embedding::
1937 sage: set_random_seed()
1938 sage: F = QuadraticField(-1, 'i')
1939 sage: M = random_matrix(F, 3)
1940 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1946 raise ValueError("the matrix 'M' must be square")
1947 if not n
.mod(2).is_zero():
1948 raise ValueError("the matrix 'M' must be a complex embedding")
1950 F
= QuadraticField(-1, 'i')
1953 # Go top-left to bottom-right (reading order), converting every
1954 # 2-by-2 block we see to a single complex element.
1956 for k
in xrange(n
/2):
1957 for j
in xrange(n
/2):
1958 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1959 if submat
[0,0] != submat
[1,1]:
1960 raise ValueError('bad on-diagonal submatrix')
1961 if submat
[0,1] != -submat
[1,0]:
1962 raise ValueError('bad off-diagonal submatrix')
1963 z
= submat
[0,0] + submat
[0,1]*i
1966 return matrix(F
, n
/2, elements
)
1969 def _embed_quaternion_matrix(M
):
1971 Embed the n-by-n quaternion matrix ``M`` into the space of real
1972 matrices of size 4n-by-4n by first sending each quaternion entry
1973 `z = a + bi + cj + dk` to the block-complex matrix
1974 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1979 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1980 sage: i,j,k = Q.gens()
1981 sage: x = 1 + 2*i + 3*j + 4*k
1982 sage: M = matrix(Q, 1, [[x]])
1983 sage: _embed_quaternion_matrix(M)
1989 Embedding is a homomorphism (isomorphism, in fact)::
1991 sage: set_random_seed()
1992 sage: n = ZZ.random_element(5)
1993 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1994 sage: X = random_matrix(Q, n)
1995 sage: Y = random_matrix(Q, n)
1996 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1997 sage: expected = _embed_quaternion_matrix(X*Y)
1998 sage: actual == expected
2002 quaternions
= M
.base_ring()
2005 raise ValueError("the matrix 'M' must be square")
2007 F
= QuadraticField(-1, 'i')
2012 t
= z
.coefficient_tuple()
2017 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2018 [-c
+ d
*i
, a
- b
*i
]])
2019 blocks
.append(_embed_complex_matrix(cplx_matrix
))
2021 # We should have real entries by now, so use the realest field
2022 # we've got for the return value.
2023 return block_matrix(quaternions
.base_ring(), n
, blocks
)
2026 def _unembed_quaternion_matrix(M
):
2028 The inverse of _embed_quaternion_matrix().
2032 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2033 ....: [-2, 1, -4, 3],
2034 ....: [-3, 4, 1, -2],
2035 ....: [-4, -3, 2, 1]])
2036 sage: _unembed_quaternion_matrix(M)
2037 [1 + 2*i + 3*j + 4*k]
2041 Unembedding is the inverse of embedding::
2043 sage: set_random_seed()
2044 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2045 sage: M = random_matrix(Q, 3)
2046 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
2052 raise ValueError("the matrix 'M' must be square")
2053 if not n
.mod(4).is_zero():
2054 raise ValueError("the matrix 'M' must be a complex embedding")
2056 Q
= QuaternionAlgebra(QQ
,-1,-1)
2059 # Go top-left to bottom-right (reading order), converting every
2060 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2063 for l
in xrange(n
/4):
2064 for m
in xrange(n
/4):
2065 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
2066 if submat
[0,0] != submat
[1,1].conjugate():
2067 raise ValueError('bad on-diagonal submatrix')
2068 if submat
[0,1] != -submat
[1,0].conjugate():
2069 raise ValueError('bad off-diagonal submatrix')
2070 z
= submat
[0,0].real() + submat
[0,0].imag()*i
2071 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
2074 return matrix(Q
, n
/4, elements
)
2077 # The usual inner product on R^n.
2079 return x
.vector().inner_product(y
.vector())
2081 # The inner product used for the real symmetric simple EJA.
2082 # We keep it as a separate function because e.g. the complex
2083 # algebra uses the same inner product, except divided by 2.
2084 def _matrix_ip(X
,Y
):
2085 X_mat
= X
.natural_representation()
2086 Y_mat
= Y
.natural_representation()
2087 return (X_mat
*Y_mat
).trace()
2090 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2092 The rank-n simple EJA consisting of real symmetric n-by-n
2093 matrices, the usual symmetric Jordan product, and the trace inner
2094 product. It has dimension `(n^2 + n)/2` over the reals.
2098 sage: J = RealSymmetricEJA(2)
2099 sage: e0, e1, e2 = J.gens()
2109 The degree of this algebra is `(n^2 + n) / 2`::
2111 sage: set_random_seed()
2112 sage: n = ZZ.random_element(1,5)
2113 sage: J = RealSymmetricEJA(n)
2114 sage: J.degree() == (n^2 + n)/2
2117 The Jordan multiplication is what we think it is::
2119 sage: set_random_seed()
2120 sage: n = ZZ.random_element(1,5)
2121 sage: J = RealSymmetricEJA(n)
2122 sage: x = J.random_element()
2123 sage: y = J.random_element()
2124 sage: actual = (x*y).natural_representation()
2125 sage: X = x.natural_representation()
2126 sage: Y = y.natural_representation()
2127 sage: expected = (X*Y + Y*X)/2
2128 sage: actual == expected
2130 sage: J(expected) == x*y
2135 def __classcall_private__(cls
, n
, field
=QQ
):
2136 S
= _real_symmetric_basis(n
, field
=field
)
2137 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2139 fdeja
= super(RealSymmetricEJA
, cls
)
2140 return fdeja
.__classcall
_private
__(cls
,
2146 def inner_product(self
, x
, y
):
2147 return _matrix_ip(x
,y
)
2150 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2152 The rank-n simple EJA consisting of complex Hermitian n-by-n
2153 matrices over the real numbers, the usual symmetric Jordan product,
2154 and the real-part-of-trace inner product. It has dimension `n^2` over
2159 The degree of this algebra is `n^2`::
2161 sage: set_random_seed()
2162 sage: n = ZZ.random_element(1,5)
2163 sage: J = ComplexHermitianEJA(n)
2164 sage: J.degree() == n^2
2167 The Jordan multiplication is what we think it is::
2169 sage: set_random_seed()
2170 sage: n = ZZ.random_element(1,5)
2171 sage: J = ComplexHermitianEJA(n)
2172 sage: x = J.random_element()
2173 sage: y = J.random_element()
2174 sage: actual = (x*y).natural_representation()
2175 sage: X = x.natural_representation()
2176 sage: Y = y.natural_representation()
2177 sage: expected = (X*Y + Y*X)/2
2178 sage: actual == expected
2180 sage: J(expected) == x*y
2185 def __classcall_private__(cls
, n
, field
=QQ
):
2186 S
= _complex_hermitian_basis(n
)
2187 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2189 fdeja
= super(ComplexHermitianEJA
, cls
)
2190 return fdeja
.__classcall
_private
__(cls
,
2196 def inner_product(self
, x
, y
):
2197 # Since a+bi on the diagonal is represented as
2202 # we'll double-count the "a" entries if we take the trace of
2204 return _matrix_ip(x
,y
)/2
2207 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2209 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2210 matrices, the usual symmetric Jordan product, and the
2211 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2216 The degree of this algebra is `n^2`::
2218 sage: set_random_seed()
2219 sage: n = ZZ.random_element(1,5)
2220 sage: J = QuaternionHermitianEJA(n)
2221 sage: J.degree() == 2*(n^2) - n
2224 The Jordan multiplication is what we think it is::
2226 sage: set_random_seed()
2227 sage: n = ZZ.random_element(1,5)
2228 sage: J = QuaternionHermitianEJA(n)
2229 sage: x = J.random_element()
2230 sage: y = J.random_element()
2231 sage: actual = (x*y).natural_representation()
2232 sage: X = x.natural_representation()
2233 sage: Y = y.natural_representation()
2234 sage: expected = (X*Y + Y*X)/2
2235 sage: actual == expected
2237 sage: J(expected) == x*y
2242 def __classcall_private__(cls
, n
, field
=QQ
):
2243 S
= _quaternion_hermitian_basis(n
)
2244 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2246 fdeja
= super(QuaternionHermitianEJA
, cls
)
2247 return fdeja
.__classcall
_private
__(cls
,
2253 def inner_product(self
, x
, y
):
2254 # Since a+bi+cj+dk on the diagonal is represented as
2256 # a + bi +cj + dk = [ a b c d]
2261 # we'll quadruple-count the "a" entries if we take the trace of
2263 return _matrix_ip(x
,y
)/4
2266 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2268 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2269 with the usual inner product and jordan product ``x*y =
2270 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2275 This multiplication table can be verified by hand::
2277 sage: J = JordanSpinEJA(4)
2278 sage: e0,e1,e2,e3 = J.gens()
2296 def __classcall_private__(cls
, n
, field
=QQ
):
2298 id_matrix
= identity_matrix(field
, n
)
2300 ei
= id_matrix
.column(i
)
2301 Qi
= zero_matrix(field
, n
)
2303 Qi
.set_column(0, ei
)
2304 Qi
+= diagonal_matrix(n
, [ei
[0]]*n
)
2305 # The addition of the diagonal matrix adds an extra ei[0] in the
2306 # upper-left corner of the matrix.
2307 Qi
[0,0] = Qi
[0,0] * ~
field(2)
2310 # The rank of the spin algebra is two, unless we're in a
2311 # one-dimensional ambient space (because the rank is bounded by
2312 # the ambient dimension).
2313 fdeja
= super(JordanSpinEJA
, cls
)
2314 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=min(n
,2))
2316 def inner_product(self
, x
, y
):
2317 return _usual_ip(x
,y
)