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
.modules
.vector_space_morphism
import VectorSpaceMorphism
11 from sage
.structure
.element
import is_Matrix
12 from sage
.structure
.category_object
import normalize_names
14 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra
import FiniteDimensionalAlgebra
15 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra_element
import FiniteDimensionalAlgebraElement
18 class FiniteDimensionalEuclideanJordanAlgebraOperator(VectorSpaceMorphism
):
19 def __init__(self
, domain_eja
, codomain_eja
, mat
):
20 # We save these so that we can output them as part of our
21 # text representation. Overriding the domain/codomain methods
22 # doesn't work because the EJAs aren't (directly) vector spaces.
23 self
._domain
_eja
= domain_eja
24 self
._codomain
_eja
= codomain_eja
26 # Otherwise, we just feed everything to the vector space morphism
28 V
= domain_eja
.vector_space()
29 W
= codomain_eja
.vector_space()
31 VectorSpaceMorphism
.__init
__(self
, homspace
, mat
)
34 def __call__(self
, x
):
36 Allow this operator to be called only on elements of an EJA.
40 sage: J = JordanSpinEJA(3)
41 sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
42 sage: id = identity_matrix(J.base_ring(), J.dimension())
43 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
48 # Overriding the single-underscore _call_ didn't work?
49 if x
not in self
._domain
_eja
:
50 raise ValueError("argument does not live in the operator's domain")
51 return self
._codomain
_eja
(self
.matrix()*x
.vector())
57 A text representation of this linear operator on a Euclidean
62 sage: J = JordanSpinEJA(2)
63 sage: id = identity_matrix(J.base_ring(), J.dimension())
64 sage: FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
65 Linear operator between finite-dimensional Euclidean Jordan
66 algebras represented by the matrix:
69 Domain: Euclidean Jordan algebra of degree 2 over Rational Field
70 Codomain: Euclidean Jordan algebra of degree 2 over Rational Field
73 msg
= ("Linear operator between finite-dimensional Euclidean Jordan "
74 "algebras represented by the matrix:\n",
78 return ''.join(msg
).format(self
.matrix(),
83 def __add__(self
, other
):
85 Add the ``other`` EJA operator to this one.
89 When we add two EJA operators, we get another one back::
91 sage: J = RealSymmetricEJA(2)
92 sage: id = identity_matrix(J.base_ring(), J.dimension())
93 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
94 sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
96 Linear operator between finite-dimensional Euclidean Jordan
97 algebras represented by the matrix:
101 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
102 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
104 If you try to add two identical vector space operators but on
105 different EJAs, that should blow up::
107 sage: J1 = RealSymmetricEJA(2)
108 sage: J2 = JordanSpinEJA(3)
109 sage: id = identity_matrix(QQ, 3)
110 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,J1,id)
111 sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,J2,id)
113 Traceback (most recent call last):
115 ValueError: operator (co)domains must match
118 if not (self
._domain
_eja
== other
._domain
_eja
and
119 self
._codomain
_eja
== other
._codomain
_eja
):
120 raise ValueError("operator (co)domains must match")
121 return FiniteDimensionalEuclideanJordanAlgebraOperator(
124 VectorSpaceMorphism
.__add
__(self
,other
))
127 def __invert__(self
):
129 Invert this EJA operator.
133 sage: J = RealSymmetricEJA(2)
134 sage: id = identity_matrix(J.base_ring(), J.dimension())
135 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
137 Linear operator between finite-dimensional Euclidean Jordan
138 algebras represented by the matrix:
142 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
143 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
146 return FiniteDimensionalEuclideanJordanAlgebraOperator(
149 VectorSpaceMorphism
.__invert
__(self
))
151 def __mul__(self
, other
):
153 Compose this EJA operator with the ``other`` one, or scale it by
154 an element of its base ring.
156 if other
in self
._codomain
_eja
.base_ring():
157 return FiniteDimensionalEuclideanJordanAlgebraOperator(
162 if not (self
._domain
_eja
== other
._codomain
_eja
):
163 raise ValueError("operator (co)domains must be compatible")
165 return FiniteDimensionalEuclideanJordanAlgebraOperator(
168 VectorSpaceMorphism
.__mul
__(self
,other
))
173 Negate this EJA operator.
177 sage: J = RealSymmetricEJA(2)
178 sage: id = identity_matrix(J.base_ring(), J.dimension())
179 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
181 Linear operator between finite-dimensional Euclidean Jordan
182 algebras represented by the matrix:
186 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
187 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
190 return FiniteDimensionalEuclideanJordanAlgebraOperator(
193 VectorSpaceMorphism
.__neg
__(self
))
196 def __pow__(self
, n
):
198 Raise this EJA operator to the power ``n``.
202 Ensure that we get back another EJA operator that can be added,
203 subtracted, et cetera::
205 sage: J = RealSymmetricEJA(2)
206 sage: id = identity_matrix(J.base_ring(), J.dimension())
207 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
208 sage: f^0 + f^1 + f^2
209 Linear operator between finite-dimensional Euclidean Jordan
210 algebras represented by the matrix:
214 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
215 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
221 # Raising a vector space morphism to the zero power gives
222 # you back a special IdentityMorphism that is useless to us.
223 rows
= self
.codomain().dimension()
224 cols
= self
.domain().dimension()
225 mat
= matrix
.identity(self
.base_ring(), rows
, cols
)
227 mat
= VectorSpaceMorphism
.__pow
__(self
,n
)
229 return FiniteDimensionalEuclideanJordanAlgebraOperator(
234 def __sub__(self
, other
):
236 Subtract ``other`` from this EJA operator.
238 return (self
+ (-other
))
241 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra
):
243 def __classcall_private__(cls
,
247 assume_associative
=False,
252 mult_table
= [b
.base_extend(field
) for b
in mult_table
]
255 if not (is_Matrix(b
) and b
.dimensions() == (n
, n
)):
256 raise ValueError("input is not a multiplication table")
257 mult_table
= tuple(mult_table
)
259 cat
= FiniteDimensionalAlgebrasWithBasis(field
)
260 cat
.or_subcategory(category
)
261 if assume_associative
:
262 cat
= cat
.Associative()
264 names
= normalize_names(n
, names
)
266 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, cls
)
267 return fda
.__classcall
__(cls
,
270 assume_associative
=assume_associative
,
274 natural_basis
=natural_basis
)
281 assume_associative
=False,
288 By definition, Jordan multiplication commutes::
290 sage: set_random_seed()
291 sage: J = random_eja()
292 sage: x = J.random_element()
293 sage: y = J.random_element()
299 self
._natural
_basis
= natural_basis
300 self
._multiplication
_table
= mult_table
301 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
310 Return a string representation of ``self``.
312 fmt
= "Euclidean Jordan algebra of degree {} over {}"
313 return fmt
.format(self
.degree(), self
.base_ring())
316 def _a_regular_element(self
):
318 Guess a regular element. Needed to compute the basis for our
319 characteristic polynomial coefficients.
322 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
323 if not z
.is_regular():
324 raise ValueError("don't know a regular element")
329 def _charpoly_basis_space(self
):
331 Return the vector space spanned by the basis used in our
332 characteristic polynomial coefficients. This is used not only to
333 compute those coefficients, but also any time we need to
334 evaluate the coefficients (like when we compute the trace or
337 z
= self
._a
_regular
_element
()
338 V
= self
.vector_space()
339 V1
= V
.span_of_basis( (z
**k
).vector() for k
in range(self
.rank()) )
340 b
= (V1
.basis() + V1
.complement().basis())
341 return V
.span_of_basis(b
)
345 def _charpoly_coeff(self
, i
):
347 Return the coefficient polynomial "a_{i}" of this algebra's
348 general characteristic polynomial.
350 Having this be a separate cached method lets us compute and
351 store the trace/determinant (a_{r-1} and a_{0} respectively)
352 separate from the entire characteristic polynomial.
354 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
355 R
= A_of_x
.base_ring()
357 # Guaranteed by theory
360 # Danger: the in-place modification is done for performance
361 # reasons (reconstructing a matrix with huge polynomial
362 # entries is slow), but I don't know how cached_method works,
363 # so it's highly possible that we're modifying some global
364 # list variable by reference, here. In other words, you
365 # probably shouldn't call this method twice on the same
366 # algebra, at the same time, in two threads
367 Ai_orig
= A_of_x
.column(i
)
368 A_of_x
.set_column(i
,xr
)
369 numerator
= A_of_x
.det()
370 A_of_x
.set_column(i
,Ai_orig
)
372 # We're relying on the theory here to ensure that each a_i is
373 # indeed back in R, and the added negative signs are to make
374 # the whole charpoly expression sum to zero.
375 return R(-numerator
/detA
)
379 def _charpoly_matrix_system(self
):
381 Compute the matrix whose entries A_ij are polynomials in
382 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
383 corresponding to `x^r` and the determinent of the matrix A =
384 [A_ij]. In other words, all of the fixed (cachable) data needed
385 to compute the coefficients of the characteristic polynomial.
390 # Construct a new algebra over a multivariate polynomial ring...
391 names
= ['X' + str(i
) for i
in range(1,n
+1)]
392 R
= PolynomialRing(self
.base_ring(), names
)
393 J
= FiniteDimensionalEuclideanJordanAlgebra(R
,
394 self
._multiplication
_table
,
397 idmat
= identity_matrix(J
.base_ring(), n
)
399 W
= self
._charpoly
_basis
_space
()
400 W
= W
.change_ring(R
.fraction_field())
402 # Starting with the standard coordinates x = (X1,X2,...,Xn)
403 # and then converting the entries to W-coordinates allows us
404 # to pass in the standard coordinates to the charpoly and get
405 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
408 # W.coordinates(x^2) eval'd at (standard z-coords)
412 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
414 # We want the middle equivalent thing in our matrix, but use
415 # the first equivalent thing instead so that we can pass in
416 # standard coordinates.
417 x
= J(vector(R
, R
.gens()))
418 l1
= [column_matrix(W
.coordinates((x
**k
).vector())) for k
in range(r
)]
419 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
420 A_of_x
= block_matrix(R
, 1, n
, (l1
+ l2
))
421 xr
= W
.coordinates((x
**r
).vector())
422 return (A_of_x
, x
, xr
, A_of_x
.det())
426 def characteristic_polynomial(self
):
431 This implementation doesn't guarantee that the polynomial
432 denominator in the coefficients is not identically zero, so
433 theoretically it could crash. The way that this is handled
434 in e.g. Faraut and Koranyi is to use a basis that guarantees
435 the denominator is non-zero. But, doing so requires knowledge
436 of at least one regular element, and we don't even know how
437 to do that. The trade-off is that, if we use the standard basis,
438 the resulting polynomial will accept the "usual" coordinates. In
439 other words, we don't have to do a change of basis before e.g.
440 computing the trace or determinant.
444 The characteristic polynomial in the spin algebra is given in
445 Alizadeh, Example 11.11::
447 sage: J = JordanSpinEJA(3)
448 sage: p = J.characteristic_polynomial(); p
449 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
450 sage: xvec = J.one().vector()
458 # The list of coefficient polynomials a_1, a_2, ..., a_n.
459 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
461 # We go to a bit of trouble here to reorder the
462 # indeterminates, so that it's easier to evaluate the
463 # characteristic polynomial at x's coordinates and get back
464 # something in terms of t, which is what we want.
466 S
= PolynomialRing(self
.base_ring(),'t')
468 S
= PolynomialRing(S
, R
.variable_names())
471 # Note: all entries past the rth should be zero. The
472 # coefficient of the highest power (x^r) is 1, but it doesn't
473 # appear in the solution vector which contains coefficients
474 # for the other powers (to make them sum to x^r).
476 a
[r
] = 1 # corresponds to x^r
478 # When the rank is equal to the dimension, trying to
479 # assign a[r] goes out-of-bounds.
480 a
.append(1) # corresponds to x^r
482 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
485 def inner_product(self
, x
, y
):
487 The inner product associated with this Euclidean Jordan algebra.
489 Defaults to the trace inner product, but can be overridden by
490 subclasses if they are sure that the necessary properties are
495 The inner product must satisfy its axiom for this algebra to truly
496 be a Euclidean Jordan Algebra::
498 sage: set_random_seed()
499 sage: J = random_eja()
500 sage: x = J.random_element()
501 sage: y = J.random_element()
502 sage: z = J.random_element()
503 sage: (x*y).inner_product(z) == y.inner_product(x*z)
507 if (not x
in self
) or (not y
in self
):
508 raise TypeError("arguments must live in this algebra")
509 return x
.trace_inner_product(y
)
512 def natural_basis(self
):
514 Return a more-natural representation of this algebra's basis.
516 Every finite-dimensional Euclidean Jordan Algebra is a direct
517 sum of five simple algebras, four of which comprise Hermitian
518 matrices. This method returns the original "natural" basis
519 for our underlying vector space. (Typically, the natural basis
520 is used to construct the multiplication table in the first place.)
522 Note that this will always return a matrix. The standard basis
523 in `R^n` will be returned as `n`-by-`1` column matrices.
527 sage: J = RealSymmetricEJA(2)
530 sage: J.natural_basis()
538 sage: J = JordanSpinEJA(2)
541 sage: J.natural_basis()
548 if self
._natural
_basis
is None:
549 return tuple( b
.vector().column() for b
in self
.basis() )
551 return self
._natural
_basis
556 Return the rank of this EJA.
558 if self
._rank
is None:
559 raise ValueError("no rank specified at genesis")
563 def vector_space(self
):
565 Return the vector space that underlies this algebra.
569 sage: J = RealSymmetricEJA(2)
570 sage: J.vector_space()
571 Vector space of dimension 3 over Rational Field
574 return self
.zero().vector().parent().ambient_vector_space()
577 class Element(FiniteDimensionalAlgebraElement
):
579 An element of a Euclidean Jordan algebra.
584 Oh man, I should not be doing this. This hides the "disabled"
585 methods ``left_matrix`` and ``matrix`` from introspection;
586 in particular it removes them from tab-completion.
588 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
589 dir(self
.__class
__) )
592 def __init__(self
, A
, elt
=None):
596 The identity in `S^n` is converted to the identity in the EJA::
598 sage: J = RealSymmetricEJA(3)
599 sage: I = identity_matrix(QQ,3)
600 sage: J(I) == J.one()
603 This skew-symmetric matrix can't be represented in the EJA::
605 sage: J = RealSymmetricEJA(3)
606 sage: A = matrix(QQ,3, lambda i,j: i-j)
608 Traceback (most recent call last):
610 ArithmeticError: vector is not in free module
613 # Goal: if we're given a matrix, and if it lives in our
614 # parent algebra's "natural ambient space," convert it
615 # into an algebra element.
617 # The catch is, we make a recursive call after converting
618 # the given matrix into a vector that lives in the algebra.
619 # This we need to try the parent class initializer first,
620 # to avoid recursing forever if we're given something that
621 # already fits into the algebra, but also happens to live
622 # in the parent's "natural ambient space" (this happens with
625 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, elt
)
627 natural_basis
= A
.natural_basis()
628 if elt
in natural_basis
[0].matrix_space():
629 # Thanks for nothing! Matrix spaces aren't vector
630 # spaces in Sage, so we have to figure out its
631 # natural-basis coordinates ourselves.
632 V
= VectorSpace(elt
.base_ring(), elt
.nrows()**2)
633 W
= V
.span( _mat2vec(s
) for s
in natural_basis
)
634 coords
= W
.coordinates(_mat2vec(elt
))
635 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, coords
)
637 def __pow__(self
, n
):
639 Return ``self`` raised to the power ``n``.
641 Jordan algebras are always power-associative; see for
642 example Faraut and Koranyi, Proposition II.1.2 (ii).
646 We have to override this because our superclass uses row vectors
647 instead of column vectors! We, on the other hand, assume column
652 sage: set_random_seed()
653 sage: x = random_eja().random_element()
654 sage: x.operator_matrix()*x.vector() == (x^2).vector()
657 A few examples of power-associativity::
659 sage: set_random_seed()
660 sage: x = random_eja().random_element()
661 sage: x*(x*x)*(x*x) == x^5
663 sage: (x*x)*(x*x*x) == x^5
666 We also know that powers operator-commute (Koecher, Chapter
669 sage: set_random_seed()
670 sage: x = random_eja().random_element()
671 sage: m = ZZ.random_element(0,10)
672 sage: n = ZZ.random_element(0,10)
673 sage: Lxm = (x^m).operator_matrix()
674 sage: Lxn = (x^n).operator_matrix()
675 sage: Lxm*Lxn == Lxn*Lxm
685 return A( (self
.operator_matrix()**(n
-1))*self
.vector() )
688 def apply_univariate_polynomial(self
, p
):
690 Apply the univariate polynomial ``p`` to this element.
692 A priori, SageMath won't allow us to apply a univariate
693 polynomial to an element of an EJA, because we don't know
694 that EJAs are rings (they are usually not associative). Of
695 course, we know that EJAs are power-associative, so the
696 operation is ultimately kosher. This function sidesteps
697 the CAS to get the answer we want and expect.
701 sage: R = PolynomialRing(QQ, 't')
703 sage: p = t^4 - t^3 + 5*t - 2
704 sage: J = RealCartesianProductEJA(5)
705 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
710 We should always get back an element of the algebra::
712 sage: set_random_seed()
713 sage: p = PolynomialRing(QQ, 't').random_element()
714 sage: J = random_eja()
715 sage: x = J.random_element()
716 sage: x.apply_univariate_polynomial(p) in J
720 if len(p
.variables()) > 1:
721 raise ValueError("not a univariate polynomial")
724 # Convert the coeficcients to the parent's base ring,
725 # because a priori they might live in an (unnecessarily)
726 # larger ring for which P.sum() would fail below.
727 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
728 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
731 def characteristic_polynomial(self
):
733 Return the characteristic polynomial of this element.
737 The rank of `R^3` is three, and the minimal polynomial of
738 the identity element is `(t-1)` from which it follows that
739 the characteristic polynomial should be `(t-1)^3`::
741 sage: J = RealCartesianProductEJA(3)
742 sage: J.one().characteristic_polynomial()
743 t^3 - 3*t^2 + 3*t - 1
745 Likewise, the characteristic of the zero element in the
746 rank-three algebra `R^{n}` should be `t^{3}`::
748 sage: J = RealCartesianProductEJA(3)
749 sage: J.zero().characteristic_polynomial()
752 The characteristic polynomial of an element should evaluate
753 to zero on that element::
755 sage: set_random_seed()
756 sage: x = RealCartesianProductEJA(3).random_element()
757 sage: p = x.characteristic_polynomial()
758 sage: x.apply_univariate_polynomial(p)
762 p
= self
.parent().characteristic_polynomial()
763 return p(*self
.vector())
766 def inner_product(self
, other
):
768 Return the parent algebra's inner product of myself and ``other``.
772 The inner product in the Jordan spin algebra is the usual
773 inner product on `R^n` (this example only works because the
774 basis for the Jordan algebra is the standard basis in `R^n`)::
776 sage: J = JordanSpinEJA(3)
777 sage: x = vector(QQ,[1,2,3])
778 sage: y = vector(QQ,[4,5,6])
779 sage: x.inner_product(y)
781 sage: J(x).inner_product(J(y))
784 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
785 multiplication is the usual matrix multiplication in `S^n`,
786 so the inner product of the identity matrix with itself
789 sage: J = RealSymmetricEJA(3)
790 sage: J.one().inner_product(J.one())
793 Likewise, the inner product on `C^n` is `<X,Y> =
794 Re(trace(X*Y))`, where we must necessarily take the real
795 part because the product of Hermitian matrices may not be
798 sage: J = ComplexHermitianEJA(3)
799 sage: J.one().inner_product(J.one())
802 Ditto for the quaternions::
804 sage: J = QuaternionHermitianEJA(3)
805 sage: J.one().inner_product(J.one())
810 Ensure that we can always compute an inner product, and that
811 it gives us back a real number::
813 sage: set_random_seed()
814 sage: J = random_eja()
815 sage: x = J.random_element()
816 sage: y = J.random_element()
817 sage: x.inner_product(y) in RR
823 raise TypeError("'other' must live in the same algebra")
825 return P
.inner_product(self
, other
)
828 def operator_commutes_with(self
, other
):
830 Return whether or not this element operator-commutes
835 The definition of a Jordan algebra says that any element
836 operator-commutes with its square::
838 sage: set_random_seed()
839 sage: x = random_eja().random_element()
840 sage: x.operator_commutes_with(x^2)
845 Test Lemma 1 from Chapter III of Koecher::
847 sage: set_random_seed()
848 sage: J = random_eja()
849 sage: u = J.random_element()
850 sage: v = J.random_element()
851 sage: lhs = u.operator_commutes_with(u*v)
852 sage: rhs = v.operator_commutes_with(u^2)
857 if not other
in self
.parent():
858 raise TypeError("'other' must live in the same algebra")
860 A
= self
.operator_matrix()
861 B
= other
.operator_matrix()
867 Return my determinant, the product of my eigenvalues.
871 sage: J = JordanSpinEJA(2)
872 sage: e0,e1 = J.gens()
873 sage: x = sum( J.gens() )
879 sage: J = JordanSpinEJA(3)
880 sage: e0,e1,e2 = J.gens()
881 sage: x = sum( J.gens() )
887 An element is invertible if and only if its determinant is
890 sage: set_random_seed()
891 sage: x = random_eja().random_element()
892 sage: x.is_invertible() == (x.det() != 0)
898 p
= P
._charpoly
_coeff
(0)
899 # The _charpoly_coeff function already adds the factor of
900 # -1 to ensure that _charpoly_coeff(0) is really what
901 # appears in front of t^{0} in the charpoly. However,
902 # we want (-1)^r times THAT for the determinant.
903 return ((-1)**r
)*p(*self
.vector())
908 Return the Jordan-multiplicative inverse of this element.
912 We appeal to the quadratic representation as in Koecher's
913 Theorem 12 in Chapter III, Section 5.
917 The inverse in the spin factor algebra is given in Alizadeh's
920 sage: set_random_seed()
921 sage: n = ZZ.random_element(1,10)
922 sage: J = JordanSpinEJA(n)
923 sage: x = J.random_element()
924 sage: while not x.is_invertible():
925 ....: x = J.random_element()
926 sage: x_vec = x.vector()
928 sage: x_bar = x_vec[1:]
929 sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
930 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
931 sage: x_inverse = coeff*inv_vec
932 sage: x.inverse() == J(x_inverse)
937 The identity element is its own inverse::
939 sage: set_random_seed()
940 sage: J = random_eja()
941 sage: J.one().inverse() == J.one()
944 If an element has an inverse, it acts like one::
946 sage: set_random_seed()
947 sage: J = random_eja()
948 sage: x = J.random_element()
949 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
952 The inverse of the inverse is what we started with::
954 sage: set_random_seed()
955 sage: J = random_eja()
956 sage: x = J.random_element()
957 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
960 The zero element is never invertible::
962 sage: set_random_seed()
963 sage: J = random_eja().zero().inverse()
964 Traceback (most recent call last):
966 ValueError: element is not invertible
969 if not self
.is_invertible():
970 raise ValueError("element is not invertible")
972 return (~self
.quadratic_representation())(self
)
975 def is_invertible(self
):
977 Return whether or not this element is invertible.
979 We can't use the superclass method because it relies on
980 the algebra being associative.
984 The usual way to do this is to check if the determinant is
985 zero, but we need the characteristic polynomial for the
986 determinant. The minimal polynomial is a lot easier to get,
987 so we use Corollary 2 in Chapter V of Koecher to check
988 whether or not the paren't algebra's zero element is a root
989 of this element's minimal polynomial.
993 The identity element is always invertible::
995 sage: set_random_seed()
996 sage: J = random_eja()
997 sage: J.one().is_invertible()
1000 The zero element is never invertible::
1002 sage: set_random_seed()
1003 sage: J = random_eja()
1004 sage: J.zero().is_invertible()
1008 zero
= self
.parent().zero()
1009 p
= self
.minimal_polynomial()
1010 return not (p(zero
) == zero
)
1013 def is_nilpotent(self
):
1015 Return whether or not some power of this element is zero.
1017 The superclass method won't work unless we're in an
1018 associative algebra, and we aren't. However, we generate
1019 an assocoative subalgebra and we're nilpotent there if and
1020 only if we're nilpotent here (probably).
1024 The identity element is never nilpotent::
1026 sage: set_random_seed()
1027 sage: random_eja().one().is_nilpotent()
1030 The additive identity is always nilpotent::
1032 sage: set_random_seed()
1033 sage: random_eja().zero().is_nilpotent()
1037 # The element we're going to call "is_nilpotent()" on.
1038 # Either myself, interpreted as an element of a finite-
1039 # dimensional algebra, or an element of an associative
1043 if self
.parent().is_associative():
1044 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1046 V
= self
.span_of_powers()
1047 assoc_subalg
= self
.subalgebra_generated_by()
1048 # Mis-design warning: the basis used for span_of_powers()
1049 # and subalgebra_generated_by() must be the same, and in
1051 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1053 # Recursive call, but should work since elt lives in an
1054 # associative algebra.
1055 return elt
.is_nilpotent()
1058 def is_regular(self
):
1060 Return whether or not this is a regular element.
1064 The identity element always has degree one, but any element
1065 linearly-independent from it is regular::
1067 sage: J = JordanSpinEJA(5)
1068 sage: J.one().is_regular()
1070 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
1071 sage: for x in J.gens():
1072 ....: (J.one() + x).is_regular()
1080 return self
.degree() == self
.parent().rank()
1085 Compute the degree of this element the straightforward way
1086 according to the definition; by appending powers to a list
1087 and figuring out its dimension (that is, whether or not
1088 they're linearly dependent).
1092 sage: J = JordanSpinEJA(4)
1093 sage: J.one().degree()
1095 sage: e0,e1,e2,e3 = J.gens()
1096 sage: (e0 - e1).degree()
1099 In the spin factor algebra (of rank two), all elements that
1100 aren't multiples of the identity are regular::
1102 sage: set_random_seed()
1103 sage: n = ZZ.random_element(1,10)
1104 sage: J = JordanSpinEJA(n)
1105 sage: x = J.random_element()
1106 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
1110 return self
.span_of_powers().dimension()
1113 def left_matrix(self
):
1115 Our parent class defines ``left_matrix`` and ``matrix``
1116 methods whose names are misleading. We don't want them.
1118 raise NotImplementedError("use operator_matrix() instead")
1120 matrix
= left_matrix
1123 def minimal_polynomial(self
):
1125 Return the minimal polynomial of this element,
1126 as a function of the variable `t`.
1130 We restrict ourselves to the associative subalgebra
1131 generated by this element, and then return the minimal
1132 polynomial of this element's operator matrix (in that
1133 subalgebra). This works by Baes Proposition 2.3.16.
1137 The minimal polynomial of the identity and zero elements are
1140 sage: set_random_seed()
1141 sage: J = random_eja()
1142 sage: J.one().minimal_polynomial()
1144 sage: J.zero().minimal_polynomial()
1147 The degree of an element is (by one definition) the degree
1148 of its minimal polynomial::
1150 sage: set_random_seed()
1151 sage: x = random_eja().random_element()
1152 sage: x.degree() == x.minimal_polynomial().degree()
1155 The minimal polynomial and the characteristic polynomial coincide
1156 and are known (see Alizadeh, Example 11.11) for all elements of
1157 the spin factor algebra that aren't scalar multiples of the
1160 sage: set_random_seed()
1161 sage: n = ZZ.random_element(2,10)
1162 sage: J = JordanSpinEJA(n)
1163 sage: y = J.random_element()
1164 sage: while y == y.coefficient(0)*J.one():
1165 ....: y = J.random_element()
1166 sage: y0 = y.vector()[0]
1167 sage: y_bar = y.vector()[1:]
1168 sage: actual = y.minimal_polynomial()
1169 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1170 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1171 sage: bool(actual == expected)
1174 The minimal polynomial should always kill its element::
1176 sage: set_random_seed()
1177 sage: x = random_eja().random_element()
1178 sage: p = x.minimal_polynomial()
1179 sage: x.apply_univariate_polynomial(p)
1183 V
= self
.span_of_powers()
1184 assoc_subalg
= self
.subalgebra_generated_by()
1185 # Mis-design warning: the basis used for span_of_powers()
1186 # and subalgebra_generated_by() must be the same, and in
1188 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1190 # We get back a symbolic polynomial in 'x' but want a real
1191 # polynomial in 't'.
1192 p_of_x
= elt
.operator_matrix().minimal_polynomial()
1193 return p_of_x
.change_variable_name('t')
1196 def natural_representation(self
):
1198 Return a more-natural representation of this element.
1200 Every finite-dimensional Euclidean Jordan Algebra is a
1201 direct sum of five simple algebras, four of which comprise
1202 Hermitian matrices. This method returns the original
1203 "natural" representation of this element as a Hermitian
1204 matrix, if it has one. If not, you get the usual representation.
1208 sage: J = ComplexHermitianEJA(3)
1211 sage: J.one().natural_representation()
1221 sage: J = QuaternionHermitianEJA(3)
1224 sage: J.one().natural_representation()
1225 [1 0 0 0 0 0 0 0 0 0 0 0]
1226 [0 1 0 0 0 0 0 0 0 0 0 0]
1227 [0 0 1 0 0 0 0 0 0 0 0 0]
1228 [0 0 0 1 0 0 0 0 0 0 0 0]
1229 [0 0 0 0 1 0 0 0 0 0 0 0]
1230 [0 0 0 0 0 1 0 0 0 0 0 0]
1231 [0 0 0 0 0 0 1 0 0 0 0 0]
1232 [0 0 0 0 0 0 0 1 0 0 0 0]
1233 [0 0 0 0 0 0 0 0 1 0 0 0]
1234 [0 0 0 0 0 0 0 0 0 1 0 0]
1235 [0 0 0 0 0 0 0 0 0 0 1 0]
1236 [0 0 0 0 0 0 0 0 0 0 0 1]
1239 B
= self
.parent().natural_basis()
1240 W
= B
[0].matrix_space()
1241 return W
.linear_combination(zip(self
.vector(), B
))
1246 Return the left-multiplication-by-this-element
1247 operator on the ambient algebra.
1251 sage: set_random_seed()
1252 sage: J = random_eja()
1253 sage: x = J.random_element()
1254 sage: y = J.random_element()
1255 sage: x.operator()(y) == x*y
1257 sage: y.operator()(x) == x*y
1262 return FiniteDimensionalEuclideanJordanAlgebraOperator(
1264 self
.operator_matrix() )
1268 def operator_matrix(self
):
1270 Return the matrix that represents left- (or right-)
1271 multiplication by this element in the parent algebra.
1273 We implement this ourselves to work around the fact that
1274 our parent class represents everything with row vectors.
1278 Test the first polarization identity from my notes, Koecher Chapter
1279 III, or from Baes (2.3)::
1281 sage: set_random_seed()
1282 sage: J = random_eja()
1283 sage: x = J.random_element()
1284 sage: y = J.random_element()
1285 sage: Lx = x.operator_matrix()
1286 sage: Ly = y.operator_matrix()
1287 sage: Lxx = (x*x).operator_matrix()
1288 sage: Lxy = (x*y).operator_matrix()
1289 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
1292 Test the second polarization identity from my notes or from
1295 sage: set_random_seed()
1296 sage: J = random_eja()
1297 sage: x = J.random_element()
1298 sage: y = J.random_element()
1299 sage: z = J.random_element()
1300 sage: Lx = x.operator_matrix()
1301 sage: Ly = y.operator_matrix()
1302 sage: Lz = z.operator_matrix()
1303 sage: Lzy = (z*y).operator_matrix()
1304 sage: Lxy = (x*y).operator_matrix()
1305 sage: Lxz = (x*z).operator_matrix()
1306 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
1309 Test the third polarization identity from my notes or from
1312 sage: set_random_seed()
1313 sage: J = random_eja()
1314 sage: u = J.random_element()
1315 sage: y = J.random_element()
1316 sage: z = J.random_element()
1317 sage: Lu = u.operator_matrix()
1318 sage: Ly = y.operator_matrix()
1319 sage: Lz = z.operator_matrix()
1320 sage: Lzy = (z*y).operator_matrix()
1321 sage: Luy = (u*y).operator_matrix()
1322 sage: Luz = (u*z).operator_matrix()
1323 sage: Luyz = (u*(y*z)).operator_matrix()
1324 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
1325 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
1326 sage: bool(lhs == rhs)
1329 Ensure that our operator's ``matrix`` method agrees with
1330 this implementation::
1332 sage: set_random_seed()
1333 sage: J = random_eja()
1334 sage: x = J.random_element()
1335 sage: x.operator().matrix() == x.operator_matrix()
1339 fda_elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1340 return fda_elt
.matrix().transpose()
1343 def quadratic_representation(self
, other
=None):
1345 Return the quadratic representation of this element.
1349 The explicit form in the spin factor algebra is given by
1350 Alizadeh's Example 11.12::
1352 sage: set_random_seed()
1353 sage: n = ZZ.random_element(1,10)
1354 sage: J = JordanSpinEJA(n)
1355 sage: x = J.random_element()
1356 sage: x_vec = x.vector()
1358 sage: x_bar = x_vec[1:]
1359 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1360 sage: B = 2*x0*x_bar.row()
1361 sage: C = 2*x0*x_bar.column()
1362 sage: D = identity_matrix(QQ, n-1)
1363 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1364 sage: D = D + 2*x_bar.tensor_product(x_bar)
1365 sage: Q = block_matrix(2,2,[A,B,C,D])
1366 sage: Q == x.quadratic_representation().matrix()
1369 Test all of the properties from Theorem 11.2 in Alizadeh::
1371 sage: set_random_seed()
1372 sage: J = random_eja()
1373 sage: x = J.random_element()
1374 sage: y = J.random_element()
1375 sage: Lx = x.operator()
1376 sage: Lxx = (x*x).operator()
1377 sage: Qx = x.quadratic_representation()
1378 sage: Qy = y.quadratic_representation()
1379 sage: Qxy = x.quadratic_representation(y)
1380 sage: Qex = J.one().quadratic_representation(x)
1381 sage: n = ZZ.random_element(10)
1382 sage: Qxn = (x^n).quadratic_representation()
1386 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1389 Property 2 (multiply on the right for :trac:`28272`):
1391 sage: alpha = QQ.random_element()
1392 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1397 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1400 sage: not x.is_invertible() or (
1403 ....: x.inverse().quadratic_representation() )
1406 sage: Qxy(J.one()) == x*y
1411 sage: not x.is_invertible() or (
1412 ....: x.quadratic_representation(x.inverse())*Qx
1413 ....: == Qx*x.quadratic_representation(x.inverse()) )
1416 sage: not x.is_invertible() or (
1417 ....: x.quadratic_representation(x.inverse())*Qx
1419 ....: 2*x.operator()*Qex - Qx )
1422 sage: 2*x.operator()*Qex - Qx == Lxx
1427 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1437 sage: not x.is_invertible() or (
1438 ....: Qx*x.inverse().operator() == Lx )
1443 sage: not x.operator_commutes_with(y) or (
1444 ....: Qx(y)^n == Qxn(y^n) )
1450 elif not other
in self
.parent():
1451 raise TypeError("'other' must live in the same algebra")
1454 M
= other
.operator()
1455 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1458 def span_of_powers(self
):
1460 Return the vector space spanned by successive powers of
1463 # The dimension of the subalgebra can't be greater than
1464 # the big algebra, so just put everything into a list
1465 # and let span() get rid of the excess.
1467 # We do the extra ambient_vector_space() in case we're messing
1468 # with polynomials and the direct parent is a module.
1469 V
= self
.parent().vector_space()
1470 return V
.span( (self
**d
).vector() for d
in xrange(V
.dimension()) )
1473 def subalgebra_generated_by(self
):
1475 Return the associative subalgebra of the parent EJA generated
1480 sage: set_random_seed()
1481 sage: x = random_eja().random_element()
1482 sage: x.subalgebra_generated_by().is_associative()
1485 Squaring in the subalgebra should be the same thing as
1486 squaring in the superalgebra::
1488 sage: set_random_seed()
1489 sage: x = random_eja().random_element()
1490 sage: u = x.subalgebra_generated_by().random_element()
1491 sage: u.operator_matrix()*u.vector() == (u**2).vector()
1495 # First get the subspace spanned by the powers of myself...
1496 V
= self
.span_of_powers()
1497 F
= self
.base_ring()
1499 # Now figure out the entries of the right-multiplication
1500 # matrix for the successive basis elements b0, b1,... of
1503 for b_right
in V
.basis():
1504 eja_b_right
= self
.parent()(b_right
)
1506 # The first row of the right-multiplication matrix by
1507 # b1 is what we get if we apply that matrix to b1. The
1508 # second row of the right multiplication matrix by b1
1509 # is what we get when we apply that matrix to b2...
1511 # IMPORTANT: this assumes that all vectors are COLUMN
1512 # vectors, unlike our superclass (which uses row vectors).
1513 for b_left
in V
.basis():
1514 eja_b_left
= self
.parent()(b_left
)
1515 # Multiply in the original EJA, but then get the
1516 # coordinates from the subalgebra in terms of its
1518 this_row
= V
.coordinates((eja_b_left
*eja_b_right
).vector())
1519 b_right_rows
.append(this_row
)
1520 b_right_matrix
= matrix(F
, b_right_rows
)
1521 mats
.append(b_right_matrix
)
1523 # It's an algebra of polynomials in one element, and EJAs
1524 # are power-associative.
1526 # TODO: choose generator names intelligently.
1527 return FiniteDimensionalEuclideanJordanAlgebra(F
, mats
, assume_associative
=True, names
='f')
1530 def subalgebra_idempotent(self
):
1532 Find an idempotent in the associative subalgebra I generate
1533 using Proposition 2.3.5 in Baes.
1537 sage: set_random_seed()
1538 sage: J = random_eja()
1539 sage: x = J.random_element()
1540 sage: while x.is_nilpotent():
1541 ....: x = J.random_element()
1542 sage: c = x.subalgebra_idempotent()
1547 if self
.is_nilpotent():
1548 raise ValueError("this only works with non-nilpotent elements!")
1550 V
= self
.span_of_powers()
1551 J
= self
.subalgebra_generated_by()
1552 # Mis-design warning: the basis used for span_of_powers()
1553 # and subalgebra_generated_by() must be the same, and in
1555 u
= J(V
.coordinates(self
.vector()))
1557 # The image of the matrix of left-u^m-multiplication
1558 # will be minimal for some natural number s...
1560 minimal_dim
= V
.dimension()
1561 for i
in xrange(1, V
.dimension()):
1562 this_dim
= (u
**i
).operator_matrix().image().dimension()
1563 if this_dim
< minimal_dim
:
1564 minimal_dim
= this_dim
1567 # Now minimal_matrix should correspond to the smallest
1568 # non-zero subspace in Baes's (or really, Koecher's)
1571 # However, we need to restrict the matrix to work on the
1572 # subspace... or do we? Can't we just solve, knowing that
1573 # A(c) = u^(s+1) should have a solution in the big space,
1576 # Beware, solve_right() means that we're using COLUMN vectors.
1577 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1579 A
= u_next
.operator_matrix()
1580 c_coordinates
= A
.solve_right(u_next
.vector())
1582 # Now c_coordinates is the idempotent we want, but it's in
1583 # the coordinate system of the subalgebra.
1585 # We need the basis for J, but as elements of the parent algebra.
1587 basis
= [self
.parent(v
) for v
in V
.basis()]
1588 return self
.parent().linear_combination(zip(c_coordinates
, basis
))
1593 Return my trace, the sum of my eigenvalues.
1597 sage: J = JordanSpinEJA(3)
1598 sage: x = sum(J.gens())
1604 sage: J = RealCartesianProductEJA(5)
1605 sage: J.one().trace()
1610 The trace of an element is a real number::
1612 sage: set_random_seed()
1613 sage: J = random_eja()
1614 sage: J.random_element().trace() in J.base_ring()
1620 p
= P
._charpoly
_coeff
(r
-1)
1621 # The _charpoly_coeff function already adds the factor of
1622 # -1 to ensure that _charpoly_coeff(r-1) is really what
1623 # appears in front of t^{r-1} in the charpoly. However,
1624 # we want the negative of THAT for the trace.
1625 return -p(*self
.vector())
1628 def trace_inner_product(self
, other
):
1630 Return the trace inner product of myself and ``other``.
1634 The trace inner product is commutative::
1636 sage: set_random_seed()
1637 sage: J = random_eja()
1638 sage: x = J.random_element(); y = J.random_element()
1639 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1642 The trace inner product is bilinear::
1644 sage: set_random_seed()
1645 sage: J = random_eja()
1646 sage: x = J.random_element()
1647 sage: y = J.random_element()
1648 sage: z = J.random_element()
1649 sage: a = QQ.random_element();
1650 sage: actual = (a*(x+z)).trace_inner_product(y)
1651 sage: expected = ( a*x.trace_inner_product(y) +
1652 ....: a*z.trace_inner_product(y) )
1653 sage: actual == expected
1655 sage: actual = x.trace_inner_product(a*(y+z))
1656 sage: expected = ( a*x.trace_inner_product(y) +
1657 ....: a*x.trace_inner_product(z) )
1658 sage: actual == expected
1661 The trace inner product satisfies the compatibility
1662 condition in the definition of a Euclidean Jordan algebra::
1664 sage: set_random_seed()
1665 sage: J = random_eja()
1666 sage: x = J.random_element()
1667 sage: y = J.random_element()
1668 sage: z = J.random_element()
1669 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1673 if not other
in self
.parent():
1674 raise TypeError("'other' must live in the same algebra")
1676 return (self
*other
).trace()
1679 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1681 Return the Euclidean Jordan Algebra corresponding to the set
1682 `R^n` under the Hadamard product.
1684 Note: this is nothing more than the Cartesian product of ``n``
1685 copies of the spin algebra. Once Cartesian product algebras
1686 are implemented, this can go.
1690 This multiplication table can be verified by hand::
1692 sage: J = RealCartesianProductEJA(3)
1693 sage: e0,e1,e2 = J.gens()
1709 def __classcall_private__(cls
, n
, field
=QQ
):
1710 # The FiniteDimensionalAlgebra constructor takes a list of
1711 # matrices, the ith representing right multiplication by the ith
1712 # basis element in the vector space. So if e_1 = (1,0,0), then
1713 # right (Hadamard) multiplication of x by e_1 picks out the first
1714 # component of x; and likewise for the ith basis element e_i.
1715 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
1716 for i
in xrange(n
) ]
1718 fdeja
= super(RealCartesianProductEJA
, cls
)
1719 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=n
)
1721 def inner_product(self
, x
, y
):
1722 return _usual_ip(x
,y
)
1727 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1731 For now, we choose a random natural number ``n`` (greater than zero)
1732 and then give you back one of the following:
1734 * The cartesian product of the rational numbers ``n`` times; this is
1735 ``QQ^n`` with the Hadamard product.
1737 * The Jordan spin algebra on ``QQ^n``.
1739 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1742 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1743 in the space of ``2n``-by-``2n`` real symmetric matrices.
1745 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1746 in the space of ``4n``-by-``4n`` real symmetric matrices.
1748 Later this might be extended to return Cartesian products of the
1754 Euclidean Jordan algebra of degree...
1758 # The max_n component lets us choose different upper bounds on the
1759 # value "n" that gets passed to the constructor. This is needed
1760 # because e.g. R^{10} is reasonable to test, while the Hermitian
1761 # 10-by-10 quaternion matrices are not.
1762 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
1764 (RealSymmetricEJA
, 5),
1765 (ComplexHermitianEJA
, 4),
1766 (QuaternionHermitianEJA
, 3)])
1767 n
= ZZ
.random_element(1, max_n
)
1768 return constructor(n
, field
=QQ
)
1772 def _real_symmetric_basis(n
, field
=QQ
):
1774 Return a basis for the space of real symmetric n-by-n matrices.
1776 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1780 for j
in xrange(i
+1):
1781 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1785 # Beware, orthogonal but not normalized!
1786 Sij
= Eij
+ Eij
.transpose()
1791 def _complex_hermitian_basis(n
, field
=QQ
):
1793 Returns a basis for the space of complex Hermitian n-by-n matrices.
1797 sage: set_random_seed()
1798 sage: n = ZZ.random_element(1,5)
1799 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1803 F
= QuadraticField(-1, 'I')
1806 # This is like the symmetric case, but we need to be careful:
1808 # * We want conjugate-symmetry, not just symmetry.
1809 # * The diagonal will (as a result) be real.
1813 for j
in xrange(i
+1):
1814 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1816 Sij
= _embed_complex_matrix(Eij
)
1819 # Beware, orthogonal but not normalized! The second one
1820 # has a minus because it's conjugated.
1821 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
1823 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
1828 def _quaternion_hermitian_basis(n
, field
=QQ
):
1830 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1834 sage: set_random_seed()
1835 sage: n = ZZ.random_element(1,5)
1836 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1840 Q
= QuaternionAlgebra(QQ
,-1,-1)
1843 # This is like the symmetric case, but we need to be careful:
1845 # * We want conjugate-symmetry, not just symmetry.
1846 # * The diagonal will (as a result) be real.
1850 for j
in xrange(i
+1):
1851 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1853 Sij
= _embed_quaternion_matrix(Eij
)
1856 # Beware, orthogonal but not normalized! The second,
1857 # third, and fourth ones have a minus because they're
1859 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
1861 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
1863 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
1865 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
1871 return vector(m
.base_ring(), m
.list())
1874 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
1876 def _multiplication_table_from_matrix_basis(basis
):
1878 At least three of the five simple Euclidean Jordan algebras have the
1879 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1880 multiplication on the right is matrix multiplication. Given a basis
1881 for the underlying matrix space, this function returns a
1882 multiplication table (obtained by looping through the basis
1883 elements) for an algebra of those matrices. A reordered copy
1884 of the basis is also returned to work around the fact that
1885 the ``span()`` in this function will change the order of the basis
1886 from what we think it is, to... something else.
1888 # In S^2, for example, we nominally have four coordinates even
1889 # though the space is of dimension three only. The vector space V
1890 # is supposed to hold the entire long vector, and the subspace W
1891 # of V will be spanned by the vectors that arise from symmetric
1892 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1893 field
= basis
[0].base_ring()
1894 dimension
= basis
[0].nrows()
1896 V
= VectorSpace(field
, dimension
**2)
1897 W
= V
.span( _mat2vec(s
) for s
in basis
)
1899 # Taking the span above reorders our basis (thanks, jerk!) so we
1900 # need to put our "matrix basis" in the same order as the
1901 # (reordered) vector basis.
1902 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
1906 # Brute force the multiplication-by-s matrix by looping
1907 # through all elements of the basis and doing the computation
1908 # to find out what the corresponding row should be. BEWARE:
1909 # these multiplication tables won't be symmetric! It therefore
1910 # becomes REALLY IMPORTANT that the underlying algebra
1911 # constructor uses ROW vectors and not COLUMN vectors. That's
1912 # why we're computing rows here and not columns.
1915 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
1916 Q_rows
.append(W
.coordinates(this_row
))
1917 Q
= matrix(field
, W
.dimension(), Q_rows
)
1923 def _embed_complex_matrix(M
):
1925 Embed the n-by-n complex matrix ``M`` into the space of real
1926 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1927 bi` to the block matrix ``[[a,b],[-b,a]]``.
1931 sage: F = QuadraticField(-1,'i')
1932 sage: x1 = F(4 - 2*i)
1933 sage: x2 = F(1 + 2*i)
1936 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1937 sage: _embed_complex_matrix(M)
1946 Embedding is a homomorphism (isomorphism, in fact)::
1948 sage: set_random_seed()
1949 sage: n = ZZ.random_element(5)
1950 sage: F = QuadraticField(-1, 'i')
1951 sage: X = random_matrix(F, n)
1952 sage: Y = random_matrix(F, n)
1953 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1954 sage: expected = _embed_complex_matrix(X*Y)
1955 sage: actual == expected
1961 raise ValueError("the matrix 'M' must be square")
1962 field
= M
.base_ring()
1967 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1969 # We can drop the imaginaries here.
1970 return block_matrix(field
.base_ring(), n
, blocks
)
1973 def _unembed_complex_matrix(M
):
1975 The inverse of _embed_complex_matrix().
1979 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1980 ....: [-2, 1, -4, 3],
1981 ....: [ 9, 10, 11, 12],
1982 ....: [-10, 9, -12, 11] ])
1983 sage: _unembed_complex_matrix(A)
1985 [ 10*i + 9 12*i + 11]
1989 Unembedding is the inverse of embedding::
1991 sage: set_random_seed()
1992 sage: F = QuadraticField(-1, 'i')
1993 sage: M = random_matrix(F, 3)
1994 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
2000 raise ValueError("the matrix 'M' must be square")
2001 if not n
.mod(2).is_zero():
2002 raise ValueError("the matrix 'M' must be a complex embedding")
2004 F
= QuadraticField(-1, 'i')
2007 # Go top-left to bottom-right (reading order), converting every
2008 # 2-by-2 block we see to a single complex element.
2010 for k
in xrange(n
/2):
2011 for j
in xrange(n
/2):
2012 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
2013 if submat
[0,0] != submat
[1,1]:
2014 raise ValueError('bad on-diagonal submatrix')
2015 if submat
[0,1] != -submat
[1,0]:
2016 raise ValueError('bad off-diagonal submatrix')
2017 z
= submat
[0,0] + submat
[0,1]*i
2020 return matrix(F
, n
/2, elements
)
2023 def _embed_quaternion_matrix(M
):
2025 Embed the n-by-n quaternion matrix ``M`` into the space of real
2026 matrices of size 4n-by-4n by first sending each quaternion entry
2027 `z = a + bi + cj + dk` to the block-complex matrix
2028 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
2033 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2034 sage: i,j,k = Q.gens()
2035 sage: x = 1 + 2*i + 3*j + 4*k
2036 sage: M = matrix(Q, 1, [[x]])
2037 sage: _embed_quaternion_matrix(M)
2043 Embedding is a homomorphism (isomorphism, in fact)::
2045 sage: set_random_seed()
2046 sage: n = ZZ.random_element(5)
2047 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2048 sage: X = random_matrix(Q, n)
2049 sage: Y = random_matrix(Q, n)
2050 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
2051 sage: expected = _embed_quaternion_matrix(X*Y)
2052 sage: actual == expected
2056 quaternions
= M
.base_ring()
2059 raise ValueError("the matrix 'M' must be square")
2061 F
= QuadraticField(-1, 'i')
2066 t
= z
.coefficient_tuple()
2071 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2072 [-c
+ d
*i
, a
- b
*i
]])
2073 blocks
.append(_embed_complex_matrix(cplx_matrix
))
2075 # We should have real entries by now, so use the realest field
2076 # we've got for the return value.
2077 return block_matrix(quaternions
.base_ring(), n
, blocks
)
2080 def _unembed_quaternion_matrix(M
):
2082 The inverse of _embed_quaternion_matrix().
2086 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2087 ....: [-2, 1, -4, 3],
2088 ....: [-3, 4, 1, -2],
2089 ....: [-4, -3, 2, 1]])
2090 sage: _unembed_quaternion_matrix(M)
2091 [1 + 2*i + 3*j + 4*k]
2095 Unembedding is the inverse of embedding::
2097 sage: set_random_seed()
2098 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2099 sage: M = random_matrix(Q, 3)
2100 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
2106 raise ValueError("the matrix 'M' must be square")
2107 if not n
.mod(4).is_zero():
2108 raise ValueError("the matrix 'M' must be a complex embedding")
2110 Q
= QuaternionAlgebra(QQ
,-1,-1)
2113 # Go top-left to bottom-right (reading order), converting every
2114 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2117 for l
in xrange(n
/4):
2118 for m
in xrange(n
/4):
2119 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
2120 if submat
[0,0] != submat
[1,1].conjugate():
2121 raise ValueError('bad on-diagonal submatrix')
2122 if submat
[0,1] != -submat
[1,0].conjugate():
2123 raise ValueError('bad off-diagonal submatrix')
2124 z
= submat
[0,0].real() + submat
[0,0].imag()*i
2125 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
2128 return matrix(Q
, n
/4, elements
)
2131 # The usual inner product on R^n.
2133 return x
.vector().inner_product(y
.vector())
2135 # The inner product used for the real symmetric simple EJA.
2136 # We keep it as a separate function because e.g. the complex
2137 # algebra uses the same inner product, except divided by 2.
2138 def _matrix_ip(X
,Y
):
2139 X_mat
= X
.natural_representation()
2140 Y_mat
= Y
.natural_representation()
2141 return (X_mat
*Y_mat
).trace()
2144 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2146 The rank-n simple EJA consisting of real symmetric n-by-n
2147 matrices, the usual symmetric Jordan product, and the trace inner
2148 product. It has dimension `(n^2 + n)/2` over the reals.
2152 sage: J = RealSymmetricEJA(2)
2153 sage: e0, e1, e2 = J.gens()
2163 The degree of this algebra is `(n^2 + n) / 2`::
2165 sage: set_random_seed()
2166 sage: n = ZZ.random_element(1,5)
2167 sage: J = RealSymmetricEJA(n)
2168 sage: J.degree() == (n^2 + n)/2
2171 The Jordan multiplication is what we think it is::
2173 sage: set_random_seed()
2174 sage: n = ZZ.random_element(1,5)
2175 sage: J = RealSymmetricEJA(n)
2176 sage: x = J.random_element()
2177 sage: y = J.random_element()
2178 sage: actual = (x*y).natural_representation()
2179 sage: X = x.natural_representation()
2180 sage: Y = y.natural_representation()
2181 sage: expected = (X*Y + Y*X)/2
2182 sage: actual == expected
2184 sage: J(expected) == x*y
2189 def __classcall_private__(cls
, n
, field
=QQ
):
2190 S
= _real_symmetric_basis(n
, field
=field
)
2191 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2193 fdeja
= super(RealSymmetricEJA
, cls
)
2194 return fdeja
.__classcall
_private
__(cls
,
2200 def inner_product(self
, x
, y
):
2201 return _matrix_ip(x
,y
)
2204 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2206 The rank-n simple EJA consisting of complex Hermitian n-by-n
2207 matrices over the real numbers, the usual symmetric Jordan product,
2208 and the real-part-of-trace inner product. It has dimension `n^2` over
2213 The degree of this algebra is `n^2`::
2215 sage: set_random_seed()
2216 sage: n = ZZ.random_element(1,5)
2217 sage: J = ComplexHermitianEJA(n)
2218 sage: J.degree() == n^2
2221 The Jordan multiplication is what we think it is::
2223 sage: set_random_seed()
2224 sage: n = ZZ.random_element(1,5)
2225 sage: J = ComplexHermitianEJA(n)
2226 sage: x = J.random_element()
2227 sage: y = J.random_element()
2228 sage: actual = (x*y).natural_representation()
2229 sage: X = x.natural_representation()
2230 sage: Y = y.natural_representation()
2231 sage: expected = (X*Y + Y*X)/2
2232 sage: actual == expected
2234 sage: J(expected) == x*y
2239 def __classcall_private__(cls
, n
, field
=QQ
):
2240 S
= _complex_hermitian_basis(n
)
2241 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2243 fdeja
= super(ComplexHermitianEJA
, cls
)
2244 return fdeja
.__classcall
_private
__(cls
,
2250 def inner_product(self
, x
, y
):
2251 # Since a+bi on the diagonal is represented as
2256 # we'll double-count the "a" entries if we take the trace of
2258 return _matrix_ip(x
,y
)/2
2261 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2263 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2264 matrices, the usual symmetric Jordan product, and the
2265 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2270 The degree of this algebra is `n^2`::
2272 sage: set_random_seed()
2273 sage: n = ZZ.random_element(1,5)
2274 sage: J = QuaternionHermitianEJA(n)
2275 sage: J.degree() == 2*(n^2) - n
2278 The Jordan multiplication is what we think it is::
2280 sage: set_random_seed()
2281 sage: n = ZZ.random_element(1,5)
2282 sage: J = QuaternionHermitianEJA(n)
2283 sage: x = J.random_element()
2284 sage: y = J.random_element()
2285 sage: actual = (x*y).natural_representation()
2286 sage: X = x.natural_representation()
2287 sage: Y = y.natural_representation()
2288 sage: expected = (X*Y + Y*X)/2
2289 sage: actual == expected
2291 sage: J(expected) == x*y
2296 def __classcall_private__(cls
, n
, field
=QQ
):
2297 S
= _quaternion_hermitian_basis(n
)
2298 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2300 fdeja
= super(QuaternionHermitianEJA
, cls
)
2301 return fdeja
.__classcall
_private
__(cls
,
2307 def inner_product(self
, x
, y
):
2308 # Since a+bi+cj+dk on the diagonal is represented as
2310 # a + bi +cj + dk = [ a b c d]
2315 # we'll quadruple-count the "a" entries if we take the trace of
2317 return _matrix_ip(x
,y
)/4
2320 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2322 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2323 with the usual inner product and jordan product ``x*y =
2324 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2329 This multiplication table can be verified by hand::
2331 sage: J = JordanSpinEJA(4)
2332 sage: e0,e1,e2,e3 = J.gens()
2350 def __classcall_private__(cls
, n
, field
=QQ
):
2352 id_matrix
= identity_matrix(field
, n
)
2354 ei
= id_matrix
.column(i
)
2355 Qi
= zero_matrix(field
, n
)
2357 Qi
.set_column(0, ei
)
2358 Qi
+= diagonal_matrix(n
, [ei
[0]]*n
)
2359 # The addition of the diagonal matrix adds an extra ei[0] in the
2360 # upper-left corner of the matrix.
2361 Qi
[0,0] = Qi
[0,0] * ~
field(2)
2364 # The rank of the spin algebra is two, unless we're in a
2365 # one-dimensional ambient space (because the rank is bounded by
2366 # the ambient dimension).
2367 fdeja
= super(JordanSpinEJA
, cls
)
2368 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=min(n
,2))
2370 def inner_product(self
, x
, y
):
2371 return _usual_ip(x
,y
)