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
.map import Map
10 from sage
.structure
.element
import is_Matrix
11 from sage
.structure
.category_object
import normalize_names
13 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra
import FiniteDimensionalAlgebra
14 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra_element
import FiniteDimensionalAlgebraElement
17 class FiniteDimensionalEuclideanJordanAlgebraOperator(Map
):
18 def __init__(self
, domain_eja
, codomain_eja
, mat
):
20 isinstance(domain_eja
, FiniteDimensionalEuclideanJordanAlgebra
) and
21 isinstance(codomain_eja
, FiniteDimensionalEuclideanJordanAlgebra
) ):
22 raise ValueError('(co)domains must be finite-dimensional Euclidean '
25 F
= domain_eja
.base_ring()
26 if not (F
== codomain_eja
.base_ring()):
27 raise ValueError("domain and codomain must have the same base ring")
29 # We need to supply something here to avoid getting the
30 # default Homset of the parent FiniteDimensionalAlgebra class,
31 # which messes up e.g. equality testing. We use FreeModules(F)
32 # instead of VectorSpaces(F) because our characteristic polynomial
33 # algorithm will need to F to be a polynomial ring at some point.
34 # When F is a field, FreeModules(F) returns VectorSpaces(F) anyway.
35 parent
= Hom(domain_eja
, codomain_eja
, FreeModules(F
))
37 # The Map initializer will set our parent to a homset, which
38 # is explicitly NOT what we want, because these ain't algebra
40 super(FiniteDimensionalEuclideanJordanAlgebraOperator
,self
).__init
__(parent
)
42 # Keep a matrix around to do all of the real work. It would
43 # be nice if we could use a VectorSpaceMorphism instead, but
44 # those use row vectors that we don't want to accidentally
45 # expose to our users.
51 Allow this operator to be called only on elements of an EJA.
55 sage: J = JordanSpinEJA(3)
56 sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
57 sage: id = identity_matrix(J.base_ring(), J.dimension())
58 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
63 return self
.codomain()(self
.matrix()*x
.vector())
66 def _add_(self
, other
):
68 Add the ``other`` EJA operator to this one.
72 When we add two EJA operators, we get another one back::
74 sage: J = RealSymmetricEJA(2)
75 sage: id = identity_matrix(J.base_ring(), J.dimension())
76 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
77 sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
79 Linear operator between finite-dimensional Euclidean Jordan
80 algebras represented by the matrix:
84 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
85 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
87 If you try to add two identical vector space operators but on
88 different EJAs, that should blow up::
90 sage: J1 = RealSymmetricEJA(2)
91 sage: J2 = JordanSpinEJA(3)
92 sage: id = identity_matrix(QQ, 3)
93 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,J1,id)
94 sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,J2,id)
96 Traceback (most recent call last):
98 TypeError: unsupported operand parent(s) for +: ...
101 return FiniteDimensionalEuclideanJordanAlgebraOperator(
104 self
.matrix() + other
.matrix())
107 def _composition_(self
, other
, homset
):
109 Compose two EJA operators to get another one (and NOT a formal
110 composite object) back.
114 sage: J1 = JordanSpinEJA(3)
115 sage: J2 = RealCartesianProductEJA(2)
116 sage: J3 = RealSymmetricEJA(1)
117 sage: mat1 = matrix(QQ, [[1,2,3],
119 sage: mat2 = matrix(QQ, [[7,8]])
120 sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,
123 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,
127 Linear operator between finite-dimensional Euclidean Jordan
128 algebras represented by the matrix:
130 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
131 Codomain: Euclidean Jordan algebra of degree 1 over Rational Field
134 return FiniteDimensionalEuclideanJordanAlgebraOperator(
137 self
.matrix()*other
.matrix())
140 def __eq__(self
, other
):
141 if self
.domain() != other
.domain():
143 if self
.codomain() != other
.codomain():
145 if self
.matrix() != other
.matrix():
149 def __invert__(self
):
151 Invert this EJA operator.
155 sage: J = RealSymmetricEJA(2)
156 sage: id = identity_matrix(J.base_ring(), J.dimension())
157 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
159 Linear operator between finite-dimensional Euclidean Jordan
160 algebras represented by the matrix:
164 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
165 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
168 return FiniteDimensionalEuclideanJordanAlgebraOperator(
174 def __mul__(self
, other
):
176 Compose two EJA operators, or scale myself by an element of the
177 ambient vector space.
179 We need to override the real ``__mul__`` function to prevent the
180 coercion framework from throwing an error when it fails to convert
181 a base ring element into a morphism.
185 We can scale an operator on a rational algebra by a rational number::
187 sage: J = RealSymmetricEJA(2)
188 sage: e0,e1,e2 = J.gens()
189 sage: x = 2*e0 + 4*e1 + 16*e2
191 Linear operator between finite-dimensional Euclidean Jordan algebras
192 represented by the matrix:
196 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
197 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
198 sage: x.operator()*(1/2)
199 Linear operator between finite-dimensional Euclidean Jordan algebras
200 represented by the matrix:
204 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
205 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
208 if other
in self
.codomain().base_ring():
209 return FiniteDimensionalEuclideanJordanAlgebraOperator(
214 # This should eventually delegate to _composition_ after performing
215 # some sanity checks for us.
216 mor
= super(FiniteDimensionalEuclideanJordanAlgebraOperator
,self
)
217 return mor
.__mul
__(other
)
222 Negate this EJA operator.
226 sage: J = RealSymmetricEJA(2)
227 sage: id = identity_matrix(J.base_ring(), J.dimension())
228 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
230 Linear operator between finite-dimensional Euclidean Jordan
231 algebras represented by the matrix:
235 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
236 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
239 return FiniteDimensionalEuclideanJordanAlgebraOperator(
245 def __pow__(self
, n
):
247 Raise this EJA operator to the power ``n``.
251 Ensure that we get back another EJA operator that can be added,
252 subtracted, et cetera::
254 sage: J = RealSymmetricEJA(2)
255 sage: id = identity_matrix(J.base_ring(), J.dimension())
256 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
257 sage: f^0 + f^1 + f^2
258 Linear operator between finite-dimensional Euclidean Jordan
259 algebras represented by the matrix:
263 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
264 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
270 # Raising a vector space morphism to the zero power gives
271 # you back a special IdentityMorphism that is useless to us.
272 rows
= self
.codomain().dimension()
273 cols
= self
.domain().dimension()
274 mat
= matrix
.identity(self
.base_ring(), rows
, cols
)
276 mat
= self
.matrix()**n
278 return FiniteDimensionalEuclideanJordanAlgebraOperator(
287 A text representation of this linear operator on a Euclidean
292 sage: J = JordanSpinEJA(2)
293 sage: id = identity_matrix(J.base_ring(), J.dimension())
294 sage: FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
295 Linear operator between finite-dimensional Euclidean Jordan
296 algebras represented by the matrix:
299 Domain: Euclidean Jordan algebra of degree 2 over Rational Field
300 Codomain: Euclidean Jordan algebra of degree 2 over Rational Field
303 msg
= ("Linear operator between finite-dimensional Euclidean Jordan "
304 "algebras represented by the matrix:\n",
308 return ''.join(msg
).format(self
.matrix(),
313 def _sub_(self
, other
):
315 Subtract ``other`` from this EJA operator.
319 sage: J = RealSymmetricEJA(2)
320 sage: id = identity_matrix(J.base_ring(),J.dimension())
321 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
323 Linear operator between finite-dimensional Euclidean Jordan
324 algebras represented by the matrix:
328 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
329 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
332 return (self
+ (-other
))
337 Return the matrix representation of this operator with respect
338 to the default bases of its (co)domain.
342 sage: J = RealSymmetricEJA(2)
343 sage: mat = matrix(J.base_ring(), J.dimension(), range(9))
344 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,mat)
354 def minimal_polynomial(self
):
356 Return the minimal polynomial of this linear operator,
357 in the variable ``t``.
361 sage: J = RealSymmetricEJA(3)
362 sage: J.one().operator().minimal_polynomial()
366 # The matrix method returns a polynomial in 'x' but want one in 't'.
367 return self
.matrix().minimal_polynomial().change_variable_name('t')
370 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra
):
372 def __classcall_private__(cls
,
376 assume_associative
=False,
381 mult_table
= [b
.base_extend(field
) for b
in mult_table
]
384 if not (is_Matrix(b
) and b
.dimensions() == (n
, n
)):
385 raise ValueError("input is not a multiplication table")
386 mult_table
= tuple(mult_table
)
388 cat
= FiniteDimensionalAlgebrasWithBasis(field
)
389 cat
.or_subcategory(category
)
390 if assume_associative
:
391 cat
= cat
.Associative()
393 names
= normalize_names(n
, names
)
395 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, cls
)
396 return fda
.__classcall
__(cls
,
399 assume_associative
=assume_associative
,
403 natural_basis
=natural_basis
)
410 assume_associative
=False,
417 By definition, Jordan multiplication commutes::
419 sage: set_random_seed()
420 sage: J = random_eja()
421 sage: x = J.random_element()
422 sage: y = J.random_element()
428 self
._natural
_basis
= natural_basis
429 self
._multiplication
_table
= mult_table
430 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
439 Return a string representation of ``self``.
441 fmt
= "Euclidean Jordan algebra of degree {} over {}"
442 return fmt
.format(self
.degree(), self
.base_ring())
445 def _a_regular_element(self
):
447 Guess a regular element. Needed to compute the basis for our
448 characteristic polynomial coefficients.
451 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
452 if not z
.is_regular():
453 raise ValueError("don't know a regular element")
458 def _charpoly_basis_space(self
):
460 Return the vector space spanned by the basis used in our
461 characteristic polynomial coefficients. This is used not only to
462 compute those coefficients, but also any time we need to
463 evaluate the coefficients (like when we compute the trace or
466 z
= self
._a
_regular
_element
()
467 V
= self
.vector_space()
468 V1
= V
.span_of_basis( (z
**k
).vector() for k
in range(self
.rank()) )
469 b
= (V1
.basis() + V1
.complement().basis())
470 return V
.span_of_basis(b
)
474 def _charpoly_coeff(self
, i
):
476 Return the coefficient polynomial "a_{i}" of this algebra's
477 general characteristic polynomial.
479 Having this be a separate cached method lets us compute and
480 store the trace/determinant (a_{r-1} and a_{0} respectively)
481 separate from the entire characteristic polynomial.
483 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
484 R
= A_of_x
.base_ring()
486 # Guaranteed by theory
489 # Danger: the in-place modification is done for performance
490 # reasons (reconstructing a matrix with huge polynomial
491 # entries is slow), but I don't know how cached_method works,
492 # so it's highly possible that we're modifying some global
493 # list variable by reference, here. In other words, you
494 # probably shouldn't call this method twice on the same
495 # algebra, at the same time, in two threads
496 Ai_orig
= A_of_x
.column(i
)
497 A_of_x
.set_column(i
,xr
)
498 numerator
= A_of_x
.det()
499 A_of_x
.set_column(i
,Ai_orig
)
501 # We're relying on the theory here to ensure that each a_i is
502 # indeed back in R, and the added negative signs are to make
503 # the whole charpoly expression sum to zero.
504 return R(-numerator
/detA
)
508 def _charpoly_matrix_system(self
):
510 Compute the matrix whose entries A_ij are polynomials in
511 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
512 corresponding to `x^r` and the determinent of the matrix A =
513 [A_ij]. In other words, all of the fixed (cachable) data needed
514 to compute the coefficients of the characteristic polynomial.
519 # Construct a new algebra over a multivariate polynomial ring...
520 names
= ['X' + str(i
) for i
in range(1,n
+1)]
521 R
= PolynomialRing(self
.base_ring(), names
)
522 J
= FiniteDimensionalEuclideanJordanAlgebra(R
,
523 self
._multiplication
_table
,
526 idmat
= identity_matrix(J
.base_ring(), n
)
528 W
= self
._charpoly
_basis
_space
()
529 W
= W
.change_ring(R
.fraction_field())
531 # Starting with the standard coordinates x = (X1,X2,...,Xn)
532 # and then converting the entries to W-coordinates allows us
533 # to pass in the standard coordinates to the charpoly and get
534 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
537 # W.coordinates(x^2) eval'd at (standard z-coords)
541 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
543 # We want the middle equivalent thing in our matrix, but use
544 # the first equivalent thing instead so that we can pass in
545 # standard coordinates.
546 x
= J(vector(R
, R
.gens()))
547 l1
= [column_matrix(W
.coordinates((x
**k
).vector())) for k
in range(r
)]
548 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
549 A_of_x
= block_matrix(R
, 1, n
, (l1
+ l2
))
550 xr
= W
.coordinates((x
**r
).vector())
551 return (A_of_x
, x
, xr
, A_of_x
.det())
555 def characteristic_polynomial(self
):
560 This implementation doesn't guarantee that the polynomial
561 denominator in the coefficients is not identically zero, so
562 theoretically it could crash. The way that this is handled
563 in e.g. Faraut and Koranyi is to use a basis that guarantees
564 the denominator is non-zero. But, doing so requires knowledge
565 of at least one regular element, and we don't even know how
566 to do that. The trade-off is that, if we use the standard basis,
567 the resulting polynomial will accept the "usual" coordinates. In
568 other words, we don't have to do a change of basis before e.g.
569 computing the trace or determinant.
573 The characteristic polynomial in the spin algebra is given in
574 Alizadeh, Example 11.11::
576 sage: J = JordanSpinEJA(3)
577 sage: p = J.characteristic_polynomial(); p
578 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
579 sage: xvec = J.one().vector()
587 # The list of coefficient polynomials a_1, a_2, ..., a_n.
588 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
590 # We go to a bit of trouble here to reorder the
591 # indeterminates, so that it's easier to evaluate the
592 # characteristic polynomial at x's coordinates and get back
593 # something in terms of t, which is what we want.
595 S
= PolynomialRing(self
.base_ring(),'t')
597 S
= PolynomialRing(S
, R
.variable_names())
600 # Note: all entries past the rth should be zero. The
601 # coefficient of the highest power (x^r) is 1, but it doesn't
602 # appear in the solution vector which contains coefficients
603 # for the other powers (to make them sum to x^r).
605 a
[r
] = 1 # corresponds to x^r
607 # When the rank is equal to the dimension, trying to
608 # assign a[r] goes out-of-bounds.
609 a
.append(1) # corresponds to x^r
611 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
614 def inner_product(self
, x
, y
):
616 The inner product associated with this Euclidean Jordan algebra.
618 Defaults to the trace inner product, but can be overridden by
619 subclasses if they are sure that the necessary properties are
624 The inner product must satisfy its axiom for this algebra to truly
625 be a Euclidean Jordan Algebra::
627 sage: set_random_seed()
628 sage: J = random_eja()
629 sage: x = J.random_element()
630 sage: y = J.random_element()
631 sage: z = J.random_element()
632 sage: (x*y).inner_product(z) == y.inner_product(x*z)
636 if (not x
in self
) or (not y
in self
):
637 raise TypeError("arguments must live in this algebra")
638 return x
.trace_inner_product(y
)
641 def natural_basis(self
):
643 Return a more-natural representation of this algebra's basis.
645 Every finite-dimensional Euclidean Jordan Algebra is a direct
646 sum of five simple algebras, four of which comprise Hermitian
647 matrices. This method returns the original "natural" basis
648 for our underlying vector space. (Typically, the natural basis
649 is used to construct the multiplication table in the first place.)
651 Note that this will always return a matrix. The standard basis
652 in `R^n` will be returned as `n`-by-`1` column matrices.
656 sage: J = RealSymmetricEJA(2)
659 sage: J.natural_basis()
667 sage: J = JordanSpinEJA(2)
670 sage: J.natural_basis()
677 if self
._natural
_basis
is None:
678 return tuple( b
.vector().column() for b
in self
.basis() )
680 return self
._natural
_basis
685 Return the rank of this EJA.
687 if self
._rank
is None:
688 raise ValueError("no rank specified at genesis")
692 def vector_space(self
):
694 Return the vector space that underlies this algebra.
698 sage: J = RealSymmetricEJA(2)
699 sage: J.vector_space()
700 Vector space of dimension 3 over Rational Field
703 return self
.zero().vector().parent().ambient_vector_space()
706 class Element(FiniteDimensionalAlgebraElement
):
708 An element of a Euclidean Jordan algebra.
713 Oh man, I should not be doing this. This hides the "disabled"
714 methods ``left_matrix`` and ``matrix`` from introspection;
715 in particular it removes them from tab-completion.
717 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
718 dir(self
.__class
__) )
721 def __init__(self
, A
, elt
=None):
725 The identity in `S^n` is converted to the identity in the EJA::
727 sage: J = RealSymmetricEJA(3)
728 sage: I = identity_matrix(QQ,3)
729 sage: J(I) == J.one()
732 This skew-symmetric matrix can't be represented in the EJA::
734 sage: J = RealSymmetricEJA(3)
735 sage: A = matrix(QQ,3, lambda i,j: i-j)
737 Traceback (most recent call last):
739 ArithmeticError: vector is not in free module
742 # Goal: if we're given a matrix, and if it lives in our
743 # parent algebra's "natural ambient space," convert it
744 # into an algebra element.
746 # The catch is, we make a recursive call after converting
747 # the given matrix into a vector that lives in the algebra.
748 # This we need to try the parent class initializer first,
749 # to avoid recursing forever if we're given something that
750 # already fits into the algebra, but also happens to live
751 # in the parent's "natural ambient space" (this happens with
754 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, elt
)
756 natural_basis
= A
.natural_basis()
757 if elt
in natural_basis
[0].matrix_space():
758 # Thanks for nothing! Matrix spaces aren't vector
759 # spaces in Sage, so we have to figure out its
760 # natural-basis coordinates ourselves.
761 V
= VectorSpace(elt
.base_ring(), elt
.nrows()**2)
762 W
= V
.span( _mat2vec(s
) for s
in natural_basis
)
763 coords
= W
.coordinates(_mat2vec(elt
))
764 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, coords
)
766 def __pow__(self
, n
):
768 Return ``self`` raised to the power ``n``.
770 Jordan algebras are always power-associative; see for
771 example Faraut and Koranyi, Proposition II.1.2 (ii).
775 We have to override this because our superclass uses row vectors
776 instead of column vectors! We, on the other hand, assume column
781 sage: set_random_seed()
782 sage: x = random_eja().random_element()
783 sage: x.operator()(x) == (x^2)
786 A few examples of power-associativity::
788 sage: set_random_seed()
789 sage: x = random_eja().random_element()
790 sage: x*(x*x)*(x*x) == x^5
792 sage: (x*x)*(x*x*x) == x^5
795 We also know that powers operator-commute (Koecher, Chapter
798 sage: set_random_seed()
799 sage: x = random_eja().random_element()
800 sage: m = ZZ.random_element(0,10)
801 sage: n = ZZ.random_element(0,10)
802 sage: Lxm = (x^m).operator()
803 sage: Lxn = (x^n).operator()
804 sage: Lxm*Lxn == Lxn*Lxm
809 return self
.parent().one()
813 return (self
.operator()**(n
-1))(self
)
816 def apply_univariate_polynomial(self
, p
):
818 Apply the univariate polynomial ``p`` to this element.
820 A priori, SageMath won't allow us to apply a univariate
821 polynomial to an element of an EJA, because we don't know
822 that EJAs are rings (they are usually not associative). Of
823 course, we know that EJAs are power-associative, so the
824 operation is ultimately kosher. This function sidesteps
825 the CAS to get the answer we want and expect.
829 sage: R = PolynomialRing(QQ, 't')
831 sage: p = t^4 - t^3 + 5*t - 2
832 sage: J = RealCartesianProductEJA(5)
833 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
838 We should always get back an element of the algebra::
840 sage: set_random_seed()
841 sage: p = PolynomialRing(QQ, 't').random_element()
842 sage: J = random_eja()
843 sage: x = J.random_element()
844 sage: x.apply_univariate_polynomial(p) in J
848 if len(p
.variables()) > 1:
849 raise ValueError("not a univariate polynomial")
852 # Convert the coeficcients to the parent's base ring,
853 # because a priori they might live in an (unnecessarily)
854 # larger ring for which P.sum() would fail below.
855 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
856 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
859 def characteristic_polynomial(self
):
861 Return the characteristic polynomial of this element.
865 The rank of `R^3` is three, and the minimal polynomial of
866 the identity element is `(t-1)` from which it follows that
867 the characteristic polynomial should be `(t-1)^3`::
869 sage: J = RealCartesianProductEJA(3)
870 sage: J.one().characteristic_polynomial()
871 t^3 - 3*t^2 + 3*t - 1
873 Likewise, the characteristic of the zero element in the
874 rank-three algebra `R^{n}` should be `t^{3}`::
876 sage: J = RealCartesianProductEJA(3)
877 sage: J.zero().characteristic_polynomial()
880 The characteristic polynomial of an element should evaluate
881 to zero on that element::
883 sage: set_random_seed()
884 sage: x = RealCartesianProductEJA(3).random_element()
885 sage: p = x.characteristic_polynomial()
886 sage: x.apply_univariate_polynomial(p)
890 p
= self
.parent().characteristic_polynomial()
891 return p(*self
.vector())
894 def inner_product(self
, other
):
896 Return the parent algebra's inner product of myself and ``other``.
900 The inner product in the Jordan spin algebra is the usual
901 inner product on `R^n` (this example only works because the
902 basis for the Jordan algebra is the standard basis in `R^n`)::
904 sage: J = JordanSpinEJA(3)
905 sage: x = vector(QQ,[1,2,3])
906 sage: y = vector(QQ,[4,5,6])
907 sage: x.inner_product(y)
909 sage: J(x).inner_product(J(y))
912 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
913 multiplication is the usual matrix multiplication in `S^n`,
914 so the inner product of the identity matrix with itself
917 sage: J = RealSymmetricEJA(3)
918 sage: J.one().inner_product(J.one())
921 Likewise, the inner product on `C^n` is `<X,Y> =
922 Re(trace(X*Y))`, where we must necessarily take the real
923 part because the product of Hermitian matrices may not be
926 sage: J = ComplexHermitianEJA(3)
927 sage: J.one().inner_product(J.one())
930 Ditto for the quaternions::
932 sage: J = QuaternionHermitianEJA(3)
933 sage: J.one().inner_product(J.one())
938 Ensure that we can always compute an inner product, and that
939 it gives us back a real number::
941 sage: set_random_seed()
942 sage: J = random_eja()
943 sage: x = J.random_element()
944 sage: y = J.random_element()
945 sage: x.inner_product(y) in RR
951 raise TypeError("'other' must live in the same algebra")
953 return P
.inner_product(self
, other
)
956 def operator_commutes_with(self
, other
):
958 Return whether or not this element operator-commutes
963 The definition of a Jordan algebra says that any element
964 operator-commutes with its square::
966 sage: set_random_seed()
967 sage: x = random_eja().random_element()
968 sage: x.operator_commutes_with(x^2)
973 Test Lemma 1 from Chapter III of Koecher::
975 sage: set_random_seed()
976 sage: J = random_eja()
977 sage: u = J.random_element()
978 sage: v = J.random_element()
979 sage: lhs = u.operator_commutes_with(u*v)
980 sage: rhs = v.operator_commutes_with(u^2)
984 Test the first polarization identity from my notes, Koecher Chapter
985 III, or from Baes (2.3)::
987 sage: set_random_seed()
988 sage: J = random_eja()
989 sage: x = J.random_element()
990 sage: y = J.random_element()
991 sage: Lx = x.operator()
992 sage: Ly = y.operator()
993 sage: Lxx = (x*x).operator()
994 sage: Lxy = (x*y).operator()
995 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
998 Test the second polarization identity from my notes or from
1001 sage: set_random_seed()
1002 sage: J = random_eja()
1003 sage: x = J.random_element()
1004 sage: y = J.random_element()
1005 sage: z = J.random_element()
1006 sage: Lx = x.operator()
1007 sage: Ly = y.operator()
1008 sage: Lz = z.operator()
1009 sage: Lzy = (z*y).operator()
1010 sage: Lxy = (x*y).operator()
1011 sage: Lxz = (x*z).operator()
1012 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
1015 Test the third polarization identity from my notes or from
1018 sage: set_random_seed()
1019 sage: J = random_eja()
1020 sage: u = J.random_element()
1021 sage: y = J.random_element()
1022 sage: z = J.random_element()
1023 sage: Lu = u.operator()
1024 sage: Ly = y.operator()
1025 sage: Lz = z.operator()
1026 sage: Lzy = (z*y).operator()
1027 sage: Luy = (u*y).operator()
1028 sage: Luz = (u*z).operator()
1029 sage: Luyz = (u*(y*z)).operator()
1030 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
1031 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
1032 sage: bool(lhs == rhs)
1036 if not other
in self
.parent():
1037 raise TypeError("'other' must live in the same algebra")
1040 B
= other
.operator()
1046 Return my determinant, the product of my eigenvalues.
1050 sage: J = JordanSpinEJA(2)
1051 sage: e0,e1 = J.gens()
1052 sage: x = sum( J.gens() )
1058 sage: J = JordanSpinEJA(3)
1059 sage: e0,e1,e2 = J.gens()
1060 sage: x = sum( J.gens() )
1066 An element is invertible if and only if its determinant is
1069 sage: set_random_seed()
1070 sage: x = random_eja().random_element()
1071 sage: x.is_invertible() == (x.det() != 0)
1077 p
= P
._charpoly
_coeff
(0)
1078 # The _charpoly_coeff function already adds the factor of
1079 # -1 to ensure that _charpoly_coeff(0) is really what
1080 # appears in front of t^{0} in the charpoly. However,
1081 # we want (-1)^r times THAT for the determinant.
1082 return ((-1)**r
)*p(*self
.vector())
1087 Return the Jordan-multiplicative inverse of this element.
1091 We appeal to the quadratic representation as in Koecher's
1092 Theorem 12 in Chapter III, Section 5.
1096 The inverse in the spin factor algebra is given in Alizadeh's
1099 sage: set_random_seed()
1100 sage: n = ZZ.random_element(1,10)
1101 sage: J = JordanSpinEJA(n)
1102 sage: x = J.random_element()
1103 sage: while not x.is_invertible():
1104 ....: x = J.random_element()
1105 sage: x_vec = x.vector()
1107 sage: x_bar = x_vec[1:]
1108 sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
1109 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
1110 sage: x_inverse = coeff*inv_vec
1111 sage: x.inverse() == J(x_inverse)
1116 The identity element is its own inverse::
1118 sage: set_random_seed()
1119 sage: J = random_eja()
1120 sage: J.one().inverse() == J.one()
1123 If an element has an inverse, it acts like one::
1125 sage: set_random_seed()
1126 sage: J = random_eja()
1127 sage: x = J.random_element()
1128 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
1131 The inverse of the inverse is what we started with::
1133 sage: set_random_seed()
1134 sage: J = random_eja()
1135 sage: x = J.random_element()
1136 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
1139 The zero element is never invertible::
1141 sage: set_random_seed()
1142 sage: J = random_eja().zero().inverse()
1143 Traceback (most recent call last):
1145 ValueError: element is not invertible
1148 if not self
.is_invertible():
1149 raise ValueError("element is not invertible")
1151 return (~self
.quadratic_representation())(self
)
1154 def is_invertible(self
):
1156 Return whether or not this element is invertible.
1158 We can't use the superclass method because it relies on
1159 the algebra being associative.
1163 The usual way to do this is to check if the determinant is
1164 zero, but we need the characteristic polynomial for the
1165 determinant. The minimal polynomial is a lot easier to get,
1166 so we use Corollary 2 in Chapter V of Koecher to check
1167 whether or not the paren't algebra's zero element is a root
1168 of this element's minimal polynomial.
1172 The identity element is always invertible::
1174 sage: set_random_seed()
1175 sage: J = random_eja()
1176 sage: J.one().is_invertible()
1179 The zero element is never invertible::
1181 sage: set_random_seed()
1182 sage: J = random_eja()
1183 sage: J.zero().is_invertible()
1187 zero
= self
.parent().zero()
1188 p
= self
.minimal_polynomial()
1189 return not (p(zero
) == zero
)
1192 def is_nilpotent(self
):
1194 Return whether or not some power of this element is zero.
1196 The superclass method won't work unless we're in an
1197 associative algebra, and we aren't. However, we generate
1198 an assocoative subalgebra and we're nilpotent there if and
1199 only if we're nilpotent here (probably).
1203 The identity element is never nilpotent::
1205 sage: set_random_seed()
1206 sage: random_eja().one().is_nilpotent()
1209 The additive identity is always nilpotent::
1211 sage: set_random_seed()
1212 sage: random_eja().zero().is_nilpotent()
1216 # The element we're going to call "is_nilpotent()" on.
1217 # Either myself, interpreted as an element of a finite-
1218 # dimensional algebra, or an element of an associative
1222 if self
.parent().is_associative():
1223 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1225 V
= self
.span_of_powers()
1226 assoc_subalg
= self
.subalgebra_generated_by()
1227 # Mis-design warning: the basis used for span_of_powers()
1228 # and subalgebra_generated_by() must be the same, and in
1230 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1232 # Recursive call, but should work since elt lives in an
1233 # associative algebra.
1234 return elt
.is_nilpotent()
1237 def is_regular(self
):
1239 Return whether or not this is a regular element.
1243 The identity element always has degree one, but any element
1244 linearly-independent from it is regular::
1246 sage: J = JordanSpinEJA(5)
1247 sage: J.one().is_regular()
1249 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
1250 sage: for x in J.gens():
1251 ....: (J.one() + x).is_regular()
1259 return self
.degree() == self
.parent().rank()
1264 Compute the degree of this element the straightforward way
1265 according to the definition; by appending powers to a list
1266 and figuring out its dimension (that is, whether or not
1267 they're linearly dependent).
1271 sage: J = JordanSpinEJA(4)
1272 sage: J.one().degree()
1274 sage: e0,e1,e2,e3 = J.gens()
1275 sage: (e0 - e1).degree()
1278 In the spin factor algebra (of rank two), all elements that
1279 aren't multiples of the identity are regular::
1281 sage: set_random_seed()
1282 sage: n = ZZ.random_element(1,10)
1283 sage: J = JordanSpinEJA(n)
1284 sage: x = J.random_element()
1285 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
1289 return self
.span_of_powers().dimension()
1292 def left_matrix(self
):
1294 Our parent class defines ``left_matrix`` and ``matrix``
1295 methods whose names are misleading. We don't want them.
1297 raise NotImplementedError("use operator().matrix() instead")
1299 matrix
= left_matrix
1302 def minimal_polynomial(self
):
1304 Return the minimal polynomial of this element,
1305 as a function of the variable `t`.
1309 We restrict ourselves to the associative subalgebra
1310 generated by this element, and then return the minimal
1311 polynomial of this element's operator matrix (in that
1312 subalgebra). This works by Baes Proposition 2.3.16.
1316 The minimal polynomial of the identity and zero elements are
1319 sage: set_random_seed()
1320 sage: J = random_eja()
1321 sage: J.one().minimal_polynomial()
1323 sage: J.zero().minimal_polynomial()
1326 The degree of an element is (by one definition) the degree
1327 of its minimal polynomial::
1329 sage: set_random_seed()
1330 sage: x = random_eja().random_element()
1331 sage: x.degree() == x.minimal_polynomial().degree()
1334 The minimal polynomial and the characteristic polynomial coincide
1335 and are known (see Alizadeh, Example 11.11) for all elements of
1336 the spin factor algebra that aren't scalar multiples of the
1339 sage: set_random_seed()
1340 sage: n = ZZ.random_element(2,10)
1341 sage: J = JordanSpinEJA(n)
1342 sage: y = J.random_element()
1343 sage: while y == y.coefficient(0)*J.one():
1344 ....: y = J.random_element()
1345 sage: y0 = y.vector()[0]
1346 sage: y_bar = y.vector()[1:]
1347 sage: actual = y.minimal_polynomial()
1348 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1349 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1350 sage: bool(actual == expected)
1353 The minimal polynomial should always kill its element::
1355 sage: set_random_seed()
1356 sage: x = random_eja().random_element()
1357 sage: p = x.minimal_polynomial()
1358 sage: x.apply_univariate_polynomial(p)
1362 V
= self
.span_of_powers()
1363 assoc_subalg
= self
.subalgebra_generated_by()
1364 # Mis-design warning: the basis used for span_of_powers()
1365 # and subalgebra_generated_by() must be the same, and in
1367 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1368 return elt
.operator().minimal_polynomial()
1372 def natural_representation(self
):
1374 Return a more-natural representation of this element.
1376 Every finite-dimensional Euclidean Jordan Algebra is a
1377 direct sum of five simple algebras, four of which comprise
1378 Hermitian matrices. This method returns the original
1379 "natural" representation of this element as a Hermitian
1380 matrix, if it has one. If not, you get the usual representation.
1384 sage: J = ComplexHermitianEJA(3)
1387 sage: J.one().natural_representation()
1397 sage: J = QuaternionHermitianEJA(3)
1400 sage: J.one().natural_representation()
1401 [1 0 0 0 0 0 0 0 0 0 0 0]
1402 [0 1 0 0 0 0 0 0 0 0 0 0]
1403 [0 0 1 0 0 0 0 0 0 0 0 0]
1404 [0 0 0 1 0 0 0 0 0 0 0 0]
1405 [0 0 0 0 1 0 0 0 0 0 0 0]
1406 [0 0 0 0 0 1 0 0 0 0 0 0]
1407 [0 0 0 0 0 0 1 0 0 0 0 0]
1408 [0 0 0 0 0 0 0 1 0 0 0 0]
1409 [0 0 0 0 0 0 0 0 1 0 0 0]
1410 [0 0 0 0 0 0 0 0 0 1 0 0]
1411 [0 0 0 0 0 0 0 0 0 0 1 0]
1412 [0 0 0 0 0 0 0 0 0 0 0 1]
1415 B
= self
.parent().natural_basis()
1416 W
= B
[0].matrix_space()
1417 return W
.linear_combination(zip(self
.vector(), B
))
1422 Return the left-multiplication-by-this-element
1423 operator on the ambient algebra.
1427 sage: set_random_seed()
1428 sage: J = random_eja()
1429 sage: x = J.random_element()
1430 sage: y = J.random_element()
1431 sage: x.operator()(y) == x*y
1433 sage: y.operator()(x) == x*y
1438 fda_elt
= FiniteDimensionalAlgebraElement(P
, self
)
1439 return FiniteDimensionalEuclideanJordanAlgebraOperator(
1442 fda_elt
.matrix().transpose() )
1445 def quadratic_representation(self
, other
=None):
1447 Return the quadratic representation of this element.
1451 The explicit form in the spin factor algebra is given by
1452 Alizadeh's Example 11.12::
1454 sage: set_random_seed()
1455 sage: n = ZZ.random_element(1,10)
1456 sage: J = JordanSpinEJA(n)
1457 sage: x = J.random_element()
1458 sage: x_vec = x.vector()
1460 sage: x_bar = x_vec[1:]
1461 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1462 sage: B = 2*x0*x_bar.row()
1463 sage: C = 2*x0*x_bar.column()
1464 sage: D = identity_matrix(QQ, n-1)
1465 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1466 sage: D = D + 2*x_bar.tensor_product(x_bar)
1467 sage: Q = block_matrix(2,2,[A,B,C,D])
1468 sage: Q == x.quadratic_representation().matrix()
1471 Test all of the properties from Theorem 11.2 in Alizadeh::
1473 sage: set_random_seed()
1474 sage: J = random_eja()
1475 sage: x = J.random_element()
1476 sage: y = J.random_element()
1477 sage: Lx = x.operator()
1478 sage: Lxx = (x*x).operator()
1479 sage: Qx = x.quadratic_representation()
1480 sage: Qy = y.quadratic_representation()
1481 sage: Qxy = x.quadratic_representation(y)
1482 sage: Qex = J.one().quadratic_representation(x)
1483 sage: n = ZZ.random_element(10)
1484 sage: Qxn = (x^n).quadratic_representation()
1488 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1491 Property 2 (multiply on the right for :trac:`28272`):
1493 sage: alpha = QQ.random_element()
1494 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1499 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1502 sage: not x.is_invertible() or (
1505 ....: x.inverse().quadratic_representation() )
1508 sage: Qxy(J.one()) == x*y
1513 sage: not x.is_invertible() or (
1514 ....: x.quadratic_representation(x.inverse())*Qx
1515 ....: == Qx*x.quadratic_representation(x.inverse()) )
1518 sage: not x.is_invertible() or (
1519 ....: x.quadratic_representation(x.inverse())*Qx
1521 ....: 2*x.operator()*Qex - Qx )
1524 sage: 2*x.operator()*Qex - Qx == Lxx
1529 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1539 sage: not x.is_invertible() or (
1540 ....: Qx*x.inverse().operator() == Lx )
1545 sage: not x.operator_commutes_with(y) or (
1546 ....: Qx(y)^n == Qxn(y^n) )
1552 elif not other
in self
.parent():
1553 raise TypeError("'other' must live in the same algebra")
1556 M
= other
.operator()
1557 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1560 def span_of_powers(self
):
1562 Return the vector space spanned by successive powers of
1565 # The dimension of the subalgebra can't be greater than
1566 # the big algebra, so just put everything into a list
1567 # and let span() get rid of the excess.
1569 # We do the extra ambient_vector_space() in case we're messing
1570 # with polynomials and the direct parent is a module.
1571 V
= self
.parent().vector_space()
1572 return V
.span( (self
**d
).vector() for d
in xrange(V
.dimension()) )
1575 def subalgebra_generated_by(self
):
1577 Return the associative subalgebra of the parent EJA generated
1582 sage: set_random_seed()
1583 sage: x = random_eja().random_element()
1584 sage: x.subalgebra_generated_by().is_associative()
1587 Squaring in the subalgebra should work the same as in
1590 sage: set_random_seed()
1591 sage: x = random_eja().random_element()
1592 sage: u = x.subalgebra_generated_by().random_element()
1593 sage: u.operator()(u) == u^2
1597 # First get the subspace spanned by the powers of myself...
1598 V
= self
.span_of_powers()
1599 F
= self
.base_ring()
1601 # Now figure out the entries of the right-multiplication
1602 # matrix for the successive basis elements b0, b1,... of
1605 for b_right
in V
.basis():
1606 eja_b_right
= self
.parent()(b_right
)
1608 # The first row of the right-multiplication matrix by
1609 # b1 is what we get if we apply that matrix to b1. The
1610 # second row of the right multiplication matrix by b1
1611 # is what we get when we apply that matrix to b2...
1613 # IMPORTANT: this assumes that all vectors are COLUMN
1614 # vectors, unlike our superclass (which uses row vectors).
1615 for b_left
in V
.basis():
1616 eja_b_left
= self
.parent()(b_left
)
1617 # Multiply in the original EJA, but then get the
1618 # coordinates from the subalgebra in terms of its
1620 this_row
= V
.coordinates((eja_b_left
*eja_b_right
).vector())
1621 b_right_rows
.append(this_row
)
1622 b_right_matrix
= matrix(F
, b_right_rows
)
1623 mats
.append(b_right_matrix
)
1625 # It's an algebra of polynomials in one element, and EJAs
1626 # are power-associative.
1628 # TODO: choose generator names intelligently.
1629 return FiniteDimensionalEuclideanJordanAlgebra(F
, mats
, assume_associative
=True, names
='f')
1632 def subalgebra_idempotent(self
):
1634 Find an idempotent in the associative subalgebra I generate
1635 using Proposition 2.3.5 in Baes.
1639 sage: set_random_seed()
1640 sage: J = random_eja()
1641 sage: x = J.random_element()
1642 sage: while x.is_nilpotent():
1643 ....: x = J.random_element()
1644 sage: c = x.subalgebra_idempotent()
1649 if self
.is_nilpotent():
1650 raise ValueError("this only works with non-nilpotent elements!")
1652 V
= self
.span_of_powers()
1653 J
= self
.subalgebra_generated_by()
1654 # Mis-design warning: the basis used for span_of_powers()
1655 # and subalgebra_generated_by() must be the same, and in
1657 u
= J(V
.coordinates(self
.vector()))
1659 # The image of the matrix of left-u^m-multiplication
1660 # will be minimal for some natural number s...
1662 minimal_dim
= V
.dimension()
1663 for i
in xrange(1, V
.dimension()):
1664 this_dim
= (u
**i
).operator().matrix().image().dimension()
1665 if this_dim
< minimal_dim
:
1666 minimal_dim
= this_dim
1669 # Now minimal_matrix should correspond to the smallest
1670 # non-zero subspace in Baes's (or really, Koecher's)
1673 # However, we need to restrict the matrix to work on the
1674 # subspace... or do we? Can't we just solve, knowing that
1675 # A(c) = u^(s+1) should have a solution in the big space,
1678 # Beware, solve_right() means that we're using COLUMN vectors.
1679 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1681 A
= u_next
.operator().matrix()
1682 c_coordinates
= A
.solve_right(u_next
.vector())
1684 # Now c_coordinates is the idempotent we want, but it's in
1685 # the coordinate system of the subalgebra.
1687 # We need the basis for J, but as elements of the parent algebra.
1689 basis
= [self
.parent(v
) for v
in V
.basis()]
1690 return self
.parent().linear_combination(zip(c_coordinates
, basis
))
1695 Return my trace, the sum of my eigenvalues.
1699 sage: J = JordanSpinEJA(3)
1700 sage: x = sum(J.gens())
1706 sage: J = RealCartesianProductEJA(5)
1707 sage: J.one().trace()
1712 The trace of an element is a real number::
1714 sage: set_random_seed()
1715 sage: J = random_eja()
1716 sage: J.random_element().trace() in J.base_ring()
1722 p
= P
._charpoly
_coeff
(r
-1)
1723 # The _charpoly_coeff function already adds the factor of
1724 # -1 to ensure that _charpoly_coeff(r-1) is really what
1725 # appears in front of t^{r-1} in the charpoly. However,
1726 # we want the negative of THAT for the trace.
1727 return -p(*self
.vector())
1730 def trace_inner_product(self
, other
):
1732 Return the trace inner product of myself and ``other``.
1736 The trace inner product is commutative::
1738 sage: set_random_seed()
1739 sage: J = random_eja()
1740 sage: x = J.random_element(); y = J.random_element()
1741 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1744 The trace inner product is bilinear::
1746 sage: set_random_seed()
1747 sage: J = random_eja()
1748 sage: x = J.random_element()
1749 sage: y = J.random_element()
1750 sage: z = J.random_element()
1751 sage: a = QQ.random_element();
1752 sage: actual = (a*(x+z)).trace_inner_product(y)
1753 sage: expected = ( a*x.trace_inner_product(y) +
1754 ....: a*z.trace_inner_product(y) )
1755 sage: actual == expected
1757 sage: actual = x.trace_inner_product(a*(y+z))
1758 sage: expected = ( a*x.trace_inner_product(y) +
1759 ....: a*x.trace_inner_product(z) )
1760 sage: actual == expected
1763 The trace inner product satisfies the compatibility
1764 condition in the definition of a Euclidean Jordan algebra::
1766 sage: set_random_seed()
1767 sage: J = random_eja()
1768 sage: x = J.random_element()
1769 sage: y = J.random_element()
1770 sage: z = J.random_element()
1771 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1775 if not other
in self
.parent():
1776 raise TypeError("'other' must live in the same algebra")
1778 return (self
*other
).trace()
1781 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1783 Return the Euclidean Jordan Algebra corresponding to the set
1784 `R^n` under the Hadamard product.
1786 Note: this is nothing more than the Cartesian product of ``n``
1787 copies of the spin algebra. Once Cartesian product algebras
1788 are implemented, this can go.
1792 This multiplication table can be verified by hand::
1794 sage: J = RealCartesianProductEJA(3)
1795 sage: e0,e1,e2 = J.gens()
1811 def __classcall_private__(cls
, n
, field
=QQ
):
1812 # The FiniteDimensionalAlgebra constructor takes a list of
1813 # matrices, the ith representing right multiplication by the ith
1814 # basis element in the vector space. So if e_1 = (1,0,0), then
1815 # right (Hadamard) multiplication of x by e_1 picks out the first
1816 # component of x; and likewise for the ith basis element e_i.
1817 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
1818 for i
in xrange(n
) ]
1820 fdeja
= super(RealCartesianProductEJA
, cls
)
1821 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=n
)
1823 def inner_product(self
, x
, y
):
1824 return _usual_ip(x
,y
)
1829 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1833 For now, we choose a random natural number ``n`` (greater than zero)
1834 and then give you back one of the following:
1836 * The cartesian product of the rational numbers ``n`` times; this is
1837 ``QQ^n`` with the Hadamard product.
1839 * The Jordan spin algebra on ``QQ^n``.
1841 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1844 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1845 in the space of ``2n``-by-``2n`` real symmetric matrices.
1847 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1848 in the space of ``4n``-by-``4n`` real symmetric matrices.
1850 Later this might be extended to return Cartesian products of the
1856 Euclidean Jordan algebra of degree...
1860 # The max_n component lets us choose different upper bounds on the
1861 # value "n" that gets passed to the constructor. This is needed
1862 # because e.g. R^{10} is reasonable to test, while the Hermitian
1863 # 10-by-10 quaternion matrices are not.
1864 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
1866 (RealSymmetricEJA
, 5),
1867 (ComplexHermitianEJA
, 4),
1868 (QuaternionHermitianEJA
, 3)])
1869 n
= ZZ
.random_element(1, max_n
)
1870 return constructor(n
, field
=QQ
)
1874 def _real_symmetric_basis(n
, field
=QQ
):
1876 Return a basis for the space of real symmetric n-by-n matrices.
1878 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1882 for j
in xrange(i
+1):
1883 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1887 # Beware, orthogonal but not normalized!
1888 Sij
= Eij
+ Eij
.transpose()
1893 def _complex_hermitian_basis(n
, field
=QQ
):
1895 Returns a basis for the space of complex Hermitian n-by-n matrices.
1899 sage: set_random_seed()
1900 sage: n = ZZ.random_element(1,5)
1901 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1905 F
= QuadraticField(-1, 'I')
1908 # This is like the symmetric case, but we need to be careful:
1910 # * We want conjugate-symmetry, not just symmetry.
1911 # * The diagonal will (as a result) be real.
1915 for j
in xrange(i
+1):
1916 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1918 Sij
= _embed_complex_matrix(Eij
)
1921 # Beware, orthogonal but not normalized! The second one
1922 # has a minus because it's conjugated.
1923 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
1925 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
1930 def _quaternion_hermitian_basis(n
, field
=QQ
):
1932 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1936 sage: set_random_seed()
1937 sage: n = ZZ.random_element(1,5)
1938 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1942 Q
= QuaternionAlgebra(QQ
,-1,-1)
1945 # This is like the symmetric case, but we need to be careful:
1947 # * We want conjugate-symmetry, not just symmetry.
1948 # * The diagonal will (as a result) be real.
1952 for j
in xrange(i
+1):
1953 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1955 Sij
= _embed_quaternion_matrix(Eij
)
1958 # Beware, orthogonal but not normalized! The second,
1959 # third, and fourth ones have a minus because they're
1961 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
1963 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
1965 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
1967 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
1973 return vector(m
.base_ring(), m
.list())
1976 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
1978 def _multiplication_table_from_matrix_basis(basis
):
1980 At least three of the five simple Euclidean Jordan algebras have the
1981 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1982 multiplication on the right is matrix multiplication. Given a basis
1983 for the underlying matrix space, this function returns a
1984 multiplication table (obtained by looping through the basis
1985 elements) for an algebra of those matrices. A reordered copy
1986 of the basis is also returned to work around the fact that
1987 the ``span()`` in this function will change the order of the basis
1988 from what we think it is, to... something else.
1990 # In S^2, for example, we nominally have four coordinates even
1991 # though the space is of dimension three only. The vector space V
1992 # is supposed to hold the entire long vector, and the subspace W
1993 # of V will be spanned by the vectors that arise from symmetric
1994 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1995 field
= basis
[0].base_ring()
1996 dimension
= basis
[0].nrows()
1998 V
= VectorSpace(field
, dimension
**2)
1999 W
= V
.span( _mat2vec(s
) for s
in basis
)
2001 # Taking the span above reorders our basis (thanks, jerk!) so we
2002 # need to put our "matrix basis" in the same order as the
2003 # (reordered) vector basis.
2004 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
2008 # Brute force the multiplication-by-s matrix by looping
2009 # through all elements of the basis and doing the computation
2010 # to find out what the corresponding row should be. BEWARE:
2011 # these multiplication tables won't be symmetric! It therefore
2012 # becomes REALLY IMPORTANT that the underlying algebra
2013 # constructor uses ROW vectors and not COLUMN vectors. That's
2014 # why we're computing rows here and not columns.
2017 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
2018 Q_rows
.append(W
.coordinates(this_row
))
2019 Q
= matrix(field
, W
.dimension(), Q_rows
)
2025 def _embed_complex_matrix(M
):
2027 Embed the n-by-n complex matrix ``M`` into the space of real
2028 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2029 bi` to the block matrix ``[[a,b],[-b,a]]``.
2033 sage: F = QuadraticField(-1,'i')
2034 sage: x1 = F(4 - 2*i)
2035 sage: x2 = F(1 + 2*i)
2038 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2039 sage: _embed_complex_matrix(M)
2048 Embedding is a homomorphism (isomorphism, in fact)::
2050 sage: set_random_seed()
2051 sage: n = ZZ.random_element(5)
2052 sage: F = QuadraticField(-1, 'i')
2053 sage: X = random_matrix(F, n)
2054 sage: Y = random_matrix(F, n)
2055 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
2056 sage: expected = _embed_complex_matrix(X*Y)
2057 sage: actual == expected
2063 raise ValueError("the matrix 'M' must be square")
2064 field
= M
.base_ring()
2069 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
2071 # We can drop the imaginaries here.
2072 return block_matrix(field
.base_ring(), n
, blocks
)
2075 def _unembed_complex_matrix(M
):
2077 The inverse of _embed_complex_matrix().
2081 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2082 ....: [-2, 1, -4, 3],
2083 ....: [ 9, 10, 11, 12],
2084 ....: [-10, 9, -12, 11] ])
2085 sage: _unembed_complex_matrix(A)
2087 [ 10*i + 9 12*i + 11]
2091 Unembedding is the inverse of embedding::
2093 sage: set_random_seed()
2094 sage: F = QuadraticField(-1, 'i')
2095 sage: M = random_matrix(F, 3)
2096 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
2102 raise ValueError("the matrix 'M' must be square")
2103 if not n
.mod(2).is_zero():
2104 raise ValueError("the matrix 'M' must be a complex embedding")
2106 F
= QuadraticField(-1, 'i')
2109 # Go top-left to bottom-right (reading order), converting every
2110 # 2-by-2 block we see to a single complex element.
2112 for k
in xrange(n
/2):
2113 for j
in xrange(n
/2):
2114 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
2115 if submat
[0,0] != submat
[1,1]:
2116 raise ValueError('bad on-diagonal submatrix')
2117 if submat
[0,1] != -submat
[1,0]:
2118 raise ValueError('bad off-diagonal submatrix')
2119 z
= submat
[0,0] + submat
[0,1]*i
2122 return matrix(F
, n
/2, elements
)
2125 def _embed_quaternion_matrix(M
):
2127 Embed the n-by-n quaternion matrix ``M`` into the space of real
2128 matrices of size 4n-by-4n by first sending each quaternion entry
2129 `z = a + bi + cj + dk` to the block-complex matrix
2130 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
2135 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2136 sage: i,j,k = Q.gens()
2137 sage: x = 1 + 2*i + 3*j + 4*k
2138 sage: M = matrix(Q, 1, [[x]])
2139 sage: _embed_quaternion_matrix(M)
2145 Embedding is a homomorphism (isomorphism, in fact)::
2147 sage: set_random_seed()
2148 sage: n = ZZ.random_element(5)
2149 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2150 sage: X = random_matrix(Q, n)
2151 sage: Y = random_matrix(Q, n)
2152 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
2153 sage: expected = _embed_quaternion_matrix(X*Y)
2154 sage: actual == expected
2158 quaternions
= M
.base_ring()
2161 raise ValueError("the matrix 'M' must be square")
2163 F
= QuadraticField(-1, 'i')
2168 t
= z
.coefficient_tuple()
2173 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2174 [-c
+ d
*i
, a
- b
*i
]])
2175 blocks
.append(_embed_complex_matrix(cplx_matrix
))
2177 # We should have real entries by now, so use the realest field
2178 # we've got for the return value.
2179 return block_matrix(quaternions
.base_ring(), n
, blocks
)
2182 def _unembed_quaternion_matrix(M
):
2184 The inverse of _embed_quaternion_matrix().
2188 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2189 ....: [-2, 1, -4, 3],
2190 ....: [-3, 4, 1, -2],
2191 ....: [-4, -3, 2, 1]])
2192 sage: _unembed_quaternion_matrix(M)
2193 [1 + 2*i + 3*j + 4*k]
2197 Unembedding is the inverse of embedding::
2199 sage: set_random_seed()
2200 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2201 sage: M = random_matrix(Q, 3)
2202 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
2208 raise ValueError("the matrix 'M' must be square")
2209 if not n
.mod(4).is_zero():
2210 raise ValueError("the matrix 'M' must be a complex embedding")
2212 Q
= QuaternionAlgebra(QQ
,-1,-1)
2215 # Go top-left to bottom-right (reading order), converting every
2216 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2219 for l
in xrange(n
/4):
2220 for m
in xrange(n
/4):
2221 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
2222 if submat
[0,0] != submat
[1,1].conjugate():
2223 raise ValueError('bad on-diagonal submatrix')
2224 if submat
[0,1] != -submat
[1,0].conjugate():
2225 raise ValueError('bad off-diagonal submatrix')
2226 z
= submat
[0,0].real() + submat
[0,0].imag()*i
2227 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
2230 return matrix(Q
, n
/4, elements
)
2233 # The usual inner product on R^n.
2235 return x
.vector().inner_product(y
.vector())
2237 # The inner product used for the real symmetric simple EJA.
2238 # We keep it as a separate function because e.g. the complex
2239 # algebra uses the same inner product, except divided by 2.
2240 def _matrix_ip(X
,Y
):
2241 X_mat
= X
.natural_representation()
2242 Y_mat
= Y
.natural_representation()
2243 return (X_mat
*Y_mat
).trace()
2246 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2248 The rank-n simple EJA consisting of real symmetric n-by-n
2249 matrices, the usual symmetric Jordan product, and the trace inner
2250 product. It has dimension `(n^2 + n)/2` over the reals.
2254 sage: J = RealSymmetricEJA(2)
2255 sage: e0, e1, e2 = J.gens()
2265 The degree of this algebra is `(n^2 + n) / 2`::
2267 sage: set_random_seed()
2268 sage: n = ZZ.random_element(1,5)
2269 sage: J = RealSymmetricEJA(n)
2270 sage: J.degree() == (n^2 + n)/2
2273 The Jordan multiplication is what we think it is::
2275 sage: set_random_seed()
2276 sage: n = ZZ.random_element(1,5)
2277 sage: J = RealSymmetricEJA(n)
2278 sage: x = J.random_element()
2279 sage: y = J.random_element()
2280 sage: actual = (x*y).natural_representation()
2281 sage: X = x.natural_representation()
2282 sage: Y = y.natural_representation()
2283 sage: expected = (X*Y + Y*X)/2
2284 sage: actual == expected
2286 sage: J(expected) == x*y
2291 def __classcall_private__(cls
, n
, field
=QQ
):
2292 S
= _real_symmetric_basis(n
, field
=field
)
2293 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2295 fdeja
= super(RealSymmetricEJA
, cls
)
2296 return fdeja
.__classcall
_private
__(cls
,
2302 def inner_product(self
, x
, y
):
2303 return _matrix_ip(x
,y
)
2306 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2308 The rank-n simple EJA consisting of complex Hermitian n-by-n
2309 matrices over the real numbers, the usual symmetric Jordan product,
2310 and the real-part-of-trace inner product. It has dimension `n^2` over
2315 The degree of this algebra is `n^2`::
2317 sage: set_random_seed()
2318 sage: n = ZZ.random_element(1,5)
2319 sage: J = ComplexHermitianEJA(n)
2320 sage: J.degree() == n^2
2323 The Jordan multiplication is what we think it is::
2325 sage: set_random_seed()
2326 sage: n = ZZ.random_element(1,5)
2327 sage: J = ComplexHermitianEJA(n)
2328 sage: x = J.random_element()
2329 sage: y = J.random_element()
2330 sage: actual = (x*y).natural_representation()
2331 sage: X = x.natural_representation()
2332 sage: Y = y.natural_representation()
2333 sage: expected = (X*Y + Y*X)/2
2334 sage: actual == expected
2336 sage: J(expected) == x*y
2341 def __classcall_private__(cls
, n
, field
=QQ
):
2342 S
= _complex_hermitian_basis(n
)
2343 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2345 fdeja
= super(ComplexHermitianEJA
, cls
)
2346 return fdeja
.__classcall
_private
__(cls
,
2352 def inner_product(self
, x
, y
):
2353 # Since a+bi on the diagonal is represented as
2358 # we'll double-count the "a" entries if we take the trace of
2360 return _matrix_ip(x
,y
)/2
2363 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2365 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2366 matrices, the usual symmetric Jordan product, and the
2367 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2372 The degree of this algebra is `n^2`::
2374 sage: set_random_seed()
2375 sage: n = ZZ.random_element(1,5)
2376 sage: J = QuaternionHermitianEJA(n)
2377 sage: J.degree() == 2*(n^2) - n
2380 The Jordan multiplication is what we think it is::
2382 sage: set_random_seed()
2383 sage: n = ZZ.random_element(1,5)
2384 sage: J = QuaternionHermitianEJA(n)
2385 sage: x = J.random_element()
2386 sage: y = J.random_element()
2387 sage: actual = (x*y).natural_representation()
2388 sage: X = x.natural_representation()
2389 sage: Y = y.natural_representation()
2390 sage: expected = (X*Y + Y*X)/2
2391 sage: actual == expected
2393 sage: J(expected) == x*y
2398 def __classcall_private__(cls
, n
, field
=QQ
):
2399 S
= _quaternion_hermitian_basis(n
)
2400 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2402 fdeja
= super(QuaternionHermitianEJA
, cls
)
2403 return fdeja
.__classcall
_private
__(cls
,
2409 def inner_product(self
, x
, y
):
2410 # Since a+bi+cj+dk on the diagonal is represented as
2412 # a + bi +cj + dk = [ a b c d]
2417 # we'll quadruple-count the "a" entries if we take the trace of
2419 return _matrix_ip(x
,y
)/4
2422 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2424 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2425 with the usual inner product and jordan product ``x*y =
2426 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2431 This multiplication table can be verified by hand::
2433 sage: J = JordanSpinEJA(4)
2434 sage: e0,e1,e2,e3 = J.gens()
2452 def __classcall_private__(cls
, n
, field
=QQ
):
2454 id_matrix
= identity_matrix(field
, n
)
2456 ei
= id_matrix
.column(i
)
2457 Qi
= zero_matrix(field
, n
)
2459 Qi
.set_column(0, ei
)
2460 Qi
+= diagonal_matrix(n
, [ei
[0]]*n
)
2461 # The addition of the diagonal matrix adds an extra ei[0] in the
2462 # upper-left corner of the matrix.
2463 Qi
[0,0] = Qi
[0,0] * ~
field(2)
2466 # The rank of the spin algebra is two, unless we're in a
2467 # one-dimensional ambient space (because the rank is bounded by
2468 # the ambient dimension).
2469 fdeja
= super(JordanSpinEJA
, cls
)
2470 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=min(n
,2))
2472 def inner_product(self
, x
, y
):
2473 return _usual_ip(x
,y
)