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():
150 def __invert__(self
):
152 Invert this EJA operator.
156 sage: J = RealSymmetricEJA(2)
157 sage: id = identity_matrix(J.base_ring(), J.dimension())
158 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
160 Linear operator between finite-dimensional Euclidean Jordan
161 algebras represented by the matrix:
165 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
166 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
169 return FiniteDimensionalEuclideanJordanAlgebraOperator(
175 def __mul__(self
, other
):
177 Compose two EJA operators, or scale myself by an element of the
178 ambient vector space.
180 We need to override the real ``__mul__`` function to prevent the
181 coercion framework from throwing an error when it fails to convert
182 a base ring element into a morphism.
186 We can scale an operator on a rational algebra by a rational number::
188 sage: J = RealSymmetricEJA(2)
189 sage: e0,e1,e2 = J.gens()
190 sage: x = 2*e0 + 4*e1 + 16*e2
192 Linear operator between finite-dimensional Euclidean Jordan algebras
193 represented by the matrix:
197 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
198 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
199 sage: x.operator()*(1/2)
200 Linear operator between finite-dimensional Euclidean Jordan algebras
201 represented by the matrix:
205 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
206 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
209 if other
in self
.codomain().base_ring():
210 return FiniteDimensionalEuclideanJordanAlgebraOperator(
215 # This should eventually delegate to _composition_ after performing
216 # some sanity checks for us.
217 mor
= super(FiniteDimensionalEuclideanJordanAlgebraOperator
,self
)
218 return mor
.__mul
__(other
)
223 Negate this EJA operator.
227 sage: J = RealSymmetricEJA(2)
228 sage: id = identity_matrix(J.base_ring(), J.dimension())
229 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
231 Linear operator between finite-dimensional Euclidean Jordan
232 algebras represented by the matrix:
236 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
237 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
240 return FiniteDimensionalEuclideanJordanAlgebraOperator(
246 def __pow__(self
, n
):
248 Raise this EJA operator to the power ``n``.
252 Ensure that we get back another EJA operator that can be added,
253 subtracted, et cetera::
255 sage: J = RealSymmetricEJA(2)
256 sage: id = identity_matrix(J.base_ring(), J.dimension())
257 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
258 sage: f^0 + f^1 + f^2
259 Linear operator between finite-dimensional Euclidean Jordan
260 algebras represented by the matrix:
264 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
265 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
271 # Raising a vector space morphism to the zero power gives
272 # you back a special IdentityMorphism that is useless to us.
273 rows
= self
.codomain().dimension()
274 cols
= self
.domain().dimension()
275 mat
= matrix
.identity(self
.base_ring(), rows
, cols
)
277 mat
= self
.matrix()**n
279 return FiniteDimensionalEuclideanJordanAlgebraOperator(
288 A text representation of this linear operator on a Euclidean
293 sage: J = JordanSpinEJA(2)
294 sage: id = identity_matrix(J.base_ring(), J.dimension())
295 sage: FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
296 Linear operator between finite-dimensional Euclidean Jordan
297 algebras represented by the matrix:
300 Domain: Euclidean Jordan algebra of degree 2 over Rational Field
301 Codomain: Euclidean Jordan algebra of degree 2 over Rational Field
304 msg
= ("Linear operator between finite-dimensional Euclidean Jordan "
305 "algebras represented by the matrix:\n",
309 return ''.join(msg
).format(self
.matrix(),
314 def _sub_(self
, other
):
316 Subtract ``other`` from this EJA operator.
320 sage: J = RealSymmetricEJA(2)
321 sage: id = identity_matrix(J.base_ring(),J.dimension())
322 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
324 Linear operator between finite-dimensional Euclidean Jordan
325 algebras represented by the matrix:
329 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
330 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
333 return (self
+ (-other
))
338 Return the matrix representation of this operator with respect
339 to the default bases of its (co)domain.
343 sage: J = RealSymmetricEJA(2)
344 sage: mat = matrix(J.base_ring(), J.dimension(), range(9))
345 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,mat)
355 def minimal_polynomial(self
):
357 Return the minimal polynomial of this linear operator,
358 in the variable ``t``.
362 sage: J = RealSymmetricEJA(3)
363 sage: J.one().operator().minimal_polynomial()
367 # The matrix method returns a polynomial in 'x' but want one in 't'.
368 return self
.matrix().minimal_polynomial().change_variable_name('t')
371 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra
):
373 def __classcall_private__(cls
,
377 assume_associative
=False,
382 mult_table
= [b
.base_extend(field
) for b
in mult_table
]
385 if not (is_Matrix(b
) and b
.dimensions() == (n
, n
)):
386 raise ValueError("input is not a multiplication table")
387 mult_table
= tuple(mult_table
)
389 cat
= FiniteDimensionalAlgebrasWithBasis(field
)
390 cat
.or_subcategory(category
)
391 if assume_associative
:
392 cat
= cat
.Associative()
394 names
= normalize_names(n
, names
)
396 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, cls
)
397 return fda
.__classcall
__(cls
,
400 assume_associative
=assume_associative
,
404 natural_basis
=natural_basis
)
411 assume_associative
=False,
418 By definition, Jordan multiplication commutes::
420 sage: set_random_seed()
421 sage: J = random_eja()
422 sage: x = J.random_element()
423 sage: y = J.random_element()
429 self
._natural
_basis
= natural_basis
430 self
._multiplication
_table
= mult_table
431 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
440 Return a string representation of ``self``.
442 fmt
= "Euclidean Jordan algebra of degree {} over {}"
443 return fmt
.format(self
.degree(), self
.base_ring())
446 def _a_regular_element(self
):
448 Guess a regular element. Needed to compute the basis for our
449 characteristic polynomial coefficients.
452 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
453 if not z
.is_regular():
454 raise ValueError("don't know a regular element")
459 def _charpoly_basis_space(self
):
461 Return the vector space spanned by the basis used in our
462 characteristic polynomial coefficients. This is used not only to
463 compute those coefficients, but also any time we need to
464 evaluate the coefficients (like when we compute the trace or
467 z
= self
._a
_regular
_element
()
468 V
= self
.vector_space()
469 V1
= V
.span_of_basis( (z
**k
).vector() for k
in range(self
.rank()) )
470 b
= (V1
.basis() + V1
.complement().basis())
471 return V
.span_of_basis(b
)
475 def _charpoly_coeff(self
, i
):
477 Return the coefficient polynomial "a_{i}" of this algebra's
478 general characteristic polynomial.
480 Having this be a separate cached method lets us compute and
481 store the trace/determinant (a_{r-1} and a_{0} respectively)
482 separate from the entire characteristic polynomial.
484 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
485 R
= A_of_x
.base_ring()
487 # Guaranteed by theory
490 # Danger: the in-place modification is done for performance
491 # reasons (reconstructing a matrix with huge polynomial
492 # entries is slow), but I don't know how cached_method works,
493 # so it's highly possible that we're modifying some global
494 # list variable by reference, here. In other words, you
495 # probably shouldn't call this method twice on the same
496 # algebra, at the same time, in two threads
497 Ai_orig
= A_of_x
.column(i
)
498 A_of_x
.set_column(i
,xr
)
499 numerator
= A_of_x
.det()
500 A_of_x
.set_column(i
,Ai_orig
)
502 # We're relying on the theory here to ensure that each a_i is
503 # indeed back in R, and the added negative signs are to make
504 # the whole charpoly expression sum to zero.
505 return R(-numerator
/detA
)
509 def _charpoly_matrix_system(self
):
511 Compute the matrix whose entries A_ij are polynomials in
512 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
513 corresponding to `x^r` and the determinent of the matrix A =
514 [A_ij]. In other words, all of the fixed (cachable) data needed
515 to compute the coefficients of the characteristic polynomial.
520 # Construct a new algebra over a multivariate polynomial ring...
521 names
= ['X' + str(i
) for i
in range(1,n
+1)]
522 R
= PolynomialRing(self
.base_ring(), names
)
523 J
= FiniteDimensionalEuclideanJordanAlgebra(R
,
524 self
._multiplication
_table
,
527 idmat
= identity_matrix(J
.base_ring(), n
)
529 W
= self
._charpoly
_basis
_space
()
530 W
= W
.change_ring(R
.fraction_field())
532 # Starting with the standard coordinates x = (X1,X2,...,Xn)
533 # and then converting the entries to W-coordinates allows us
534 # to pass in the standard coordinates to the charpoly and get
535 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
538 # W.coordinates(x^2) eval'd at (standard z-coords)
542 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
544 # We want the middle equivalent thing in our matrix, but use
545 # the first equivalent thing instead so that we can pass in
546 # standard coordinates.
547 x
= J(vector(R
, R
.gens()))
548 l1
= [column_matrix(W
.coordinates((x
**k
).vector())) for k
in range(r
)]
549 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
550 A_of_x
= block_matrix(R
, 1, n
, (l1
+ l2
))
551 xr
= W
.coordinates((x
**r
).vector())
552 return (A_of_x
, x
, xr
, A_of_x
.det())
556 def characteristic_polynomial(self
):
561 This implementation doesn't guarantee that the polynomial
562 denominator in the coefficients is not identically zero, so
563 theoretically it could crash. The way that this is handled
564 in e.g. Faraut and Koranyi is to use a basis that guarantees
565 the denominator is non-zero. But, doing so requires knowledge
566 of at least one regular element, and we don't even know how
567 to do that. The trade-off is that, if we use the standard basis,
568 the resulting polynomial will accept the "usual" coordinates. In
569 other words, we don't have to do a change of basis before e.g.
570 computing the trace or determinant.
574 The characteristic polynomial in the spin algebra is given in
575 Alizadeh, Example 11.11::
577 sage: J = JordanSpinEJA(3)
578 sage: p = J.characteristic_polynomial(); p
579 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
580 sage: xvec = J.one().vector()
588 # The list of coefficient polynomials a_1, a_2, ..., a_n.
589 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
591 # We go to a bit of trouble here to reorder the
592 # indeterminates, so that it's easier to evaluate the
593 # characteristic polynomial at x's coordinates and get back
594 # something in terms of t, which is what we want.
596 S
= PolynomialRing(self
.base_ring(),'t')
598 S
= PolynomialRing(S
, R
.variable_names())
601 # Note: all entries past the rth should be zero. The
602 # coefficient of the highest power (x^r) is 1, but it doesn't
603 # appear in the solution vector which contains coefficients
604 # for the other powers (to make them sum to x^r).
606 a
[r
] = 1 # corresponds to x^r
608 # When the rank is equal to the dimension, trying to
609 # assign a[r] goes out-of-bounds.
610 a
.append(1) # corresponds to x^r
612 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
615 def inner_product(self
, x
, y
):
617 The inner product associated with this Euclidean Jordan algebra.
619 Defaults to the trace inner product, but can be overridden by
620 subclasses if they are sure that the necessary properties are
625 The inner product must satisfy its axiom for this algebra to truly
626 be a Euclidean Jordan Algebra::
628 sage: set_random_seed()
629 sage: J = random_eja()
630 sage: x = J.random_element()
631 sage: y = J.random_element()
632 sage: z = J.random_element()
633 sage: (x*y).inner_product(z) == y.inner_product(x*z)
637 if (not x
in self
) or (not y
in self
):
638 raise TypeError("arguments must live in this algebra")
639 return x
.trace_inner_product(y
)
642 def natural_basis(self
):
644 Return a more-natural representation of this algebra's basis.
646 Every finite-dimensional Euclidean Jordan Algebra is a direct
647 sum of five simple algebras, four of which comprise Hermitian
648 matrices. This method returns the original "natural" basis
649 for our underlying vector space. (Typically, the natural basis
650 is used to construct the multiplication table in the first place.)
652 Note that this will always return a matrix. The standard basis
653 in `R^n` will be returned as `n`-by-`1` column matrices.
657 sage: J = RealSymmetricEJA(2)
660 sage: J.natural_basis()
668 sage: J = JordanSpinEJA(2)
671 sage: J.natural_basis()
678 if self
._natural
_basis
is None:
679 return tuple( b
.vector().column() for b
in self
.basis() )
681 return self
._natural
_basis
686 Return the rank of this EJA.
688 if self
._rank
is None:
689 raise ValueError("no rank specified at genesis")
693 def vector_space(self
):
695 Return the vector space that underlies this algebra.
699 sage: J = RealSymmetricEJA(2)
700 sage: J.vector_space()
701 Vector space of dimension 3 over Rational Field
704 return self
.zero().vector().parent().ambient_vector_space()
707 class Element(FiniteDimensionalAlgebraElement
):
709 An element of a Euclidean Jordan algebra.
714 Oh man, I should not be doing this. This hides the "disabled"
715 methods ``left_matrix`` and ``matrix`` from introspection;
716 in particular it removes them from tab-completion.
718 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
719 dir(self
.__class
__) )
722 def __init__(self
, A
, elt
=None):
726 The identity in `S^n` is converted to the identity in the EJA::
728 sage: J = RealSymmetricEJA(3)
729 sage: I = identity_matrix(QQ,3)
730 sage: J(I) == J.one()
733 This skew-symmetric matrix can't be represented in the EJA::
735 sage: J = RealSymmetricEJA(3)
736 sage: A = matrix(QQ,3, lambda i,j: i-j)
738 Traceback (most recent call last):
740 ArithmeticError: vector is not in free module
743 # Goal: if we're given a matrix, and if it lives in our
744 # parent algebra's "natural ambient space," convert it
745 # into an algebra element.
747 # The catch is, we make a recursive call after converting
748 # the given matrix into a vector that lives in the algebra.
749 # This we need to try the parent class initializer first,
750 # to avoid recursing forever if we're given something that
751 # already fits into the algebra, but also happens to live
752 # in the parent's "natural ambient space" (this happens with
755 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, elt
)
757 natural_basis
= A
.natural_basis()
758 if elt
in natural_basis
[0].matrix_space():
759 # Thanks for nothing! Matrix spaces aren't vector
760 # spaces in Sage, so we have to figure out its
761 # natural-basis coordinates ourselves.
762 V
= VectorSpace(elt
.base_ring(), elt
.nrows()**2)
763 W
= V
.span( _mat2vec(s
) for s
in natural_basis
)
764 coords
= W
.coordinates(_mat2vec(elt
))
765 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, coords
)
767 def __pow__(self
, n
):
769 Return ``self`` raised to the power ``n``.
771 Jordan algebras are always power-associative; see for
772 example Faraut and Koranyi, Proposition II.1.2 (ii).
776 We have to override this because our superclass uses row vectors
777 instead of column vectors! We, on the other hand, assume column
782 sage: set_random_seed()
783 sage: x = random_eja().random_element()
784 sage: x.operator()(x) == (x^2)
787 A few examples of power-associativity::
789 sage: set_random_seed()
790 sage: x = random_eja().random_element()
791 sage: x*(x*x)*(x*x) == x^5
793 sage: (x*x)*(x*x*x) == x^5
796 We also know that powers operator-commute (Koecher, Chapter
799 sage: set_random_seed()
800 sage: x = random_eja().random_element()
801 sage: m = ZZ.random_element(0,10)
802 sage: n = ZZ.random_element(0,10)
803 sage: Lxm = (x^m).operator()
804 sage: Lxn = (x^n).operator()
805 sage: Lxm*Lxn == Lxn*Lxm
810 return self
.parent().one()
814 return (self
.operator()**(n
-1))(self
)
817 def apply_univariate_polynomial(self
, p
):
819 Apply the univariate polynomial ``p`` to this element.
821 A priori, SageMath won't allow us to apply a univariate
822 polynomial to an element of an EJA, because we don't know
823 that EJAs are rings (they are usually not associative). Of
824 course, we know that EJAs are power-associative, so the
825 operation is ultimately kosher. This function sidesteps
826 the CAS to get the answer we want and expect.
830 sage: R = PolynomialRing(QQ, 't')
832 sage: p = t^4 - t^3 + 5*t - 2
833 sage: J = RealCartesianProductEJA(5)
834 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
839 We should always get back an element of the algebra::
841 sage: set_random_seed()
842 sage: p = PolynomialRing(QQ, 't').random_element()
843 sage: J = random_eja()
844 sage: x = J.random_element()
845 sage: x.apply_univariate_polynomial(p) in J
849 if len(p
.variables()) > 1:
850 raise ValueError("not a univariate polynomial")
853 # Convert the coeficcients to the parent's base ring,
854 # because a priori they might live in an (unnecessarily)
855 # larger ring for which P.sum() would fail below.
856 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
857 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
860 def characteristic_polynomial(self
):
862 Return the characteristic polynomial of this element.
866 The rank of `R^3` is three, and the minimal polynomial of
867 the identity element is `(t-1)` from which it follows that
868 the characteristic polynomial should be `(t-1)^3`::
870 sage: J = RealCartesianProductEJA(3)
871 sage: J.one().characteristic_polynomial()
872 t^3 - 3*t^2 + 3*t - 1
874 Likewise, the characteristic of the zero element in the
875 rank-three algebra `R^{n}` should be `t^{3}`::
877 sage: J = RealCartesianProductEJA(3)
878 sage: J.zero().characteristic_polynomial()
881 The characteristic polynomial of an element should evaluate
882 to zero on that element::
884 sage: set_random_seed()
885 sage: x = RealCartesianProductEJA(3).random_element()
886 sage: p = x.characteristic_polynomial()
887 sage: x.apply_univariate_polynomial(p)
891 p
= self
.parent().characteristic_polynomial()
892 return p(*self
.vector())
895 def inner_product(self
, other
):
897 Return the parent algebra's inner product of myself and ``other``.
901 The inner product in the Jordan spin algebra is the usual
902 inner product on `R^n` (this example only works because the
903 basis for the Jordan algebra is the standard basis in `R^n`)::
905 sage: J = JordanSpinEJA(3)
906 sage: x = vector(QQ,[1,2,3])
907 sage: y = vector(QQ,[4,5,6])
908 sage: x.inner_product(y)
910 sage: J(x).inner_product(J(y))
913 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
914 multiplication is the usual matrix multiplication in `S^n`,
915 so the inner product of the identity matrix with itself
918 sage: J = RealSymmetricEJA(3)
919 sage: J.one().inner_product(J.one())
922 Likewise, the inner product on `C^n` is `<X,Y> =
923 Re(trace(X*Y))`, where we must necessarily take the real
924 part because the product of Hermitian matrices may not be
927 sage: J = ComplexHermitianEJA(3)
928 sage: J.one().inner_product(J.one())
931 Ditto for the quaternions::
933 sage: J = QuaternionHermitianEJA(3)
934 sage: J.one().inner_product(J.one())
939 Ensure that we can always compute an inner product, and that
940 it gives us back a real number::
942 sage: set_random_seed()
943 sage: J = random_eja()
944 sage: x = J.random_element()
945 sage: y = J.random_element()
946 sage: x.inner_product(y) in RR
952 raise TypeError("'other' must live in the same algebra")
954 return P
.inner_product(self
, other
)
957 def operator_commutes_with(self
, other
):
959 Return whether or not this element operator-commutes
964 The definition of a Jordan algebra says that any element
965 operator-commutes with its square::
967 sage: set_random_seed()
968 sage: x = random_eja().random_element()
969 sage: x.operator_commutes_with(x^2)
974 Test Lemma 1 from Chapter III of Koecher::
976 sage: set_random_seed()
977 sage: J = random_eja()
978 sage: u = J.random_element()
979 sage: v = J.random_element()
980 sage: lhs = u.operator_commutes_with(u*v)
981 sage: rhs = v.operator_commutes_with(u^2)
985 Test the first polarization identity from my notes, Koecher Chapter
986 III, or from Baes (2.3)::
988 sage: set_random_seed()
989 sage: J = random_eja()
990 sage: x = J.random_element()
991 sage: y = J.random_element()
992 sage: Lx = x.operator()
993 sage: Ly = y.operator()
994 sage: Lxx = (x*x).operator()
995 sage: Lxy = (x*y).operator()
996 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
999 Test the second polarization identity from my notes or from
1002 sage: set_random_seed()
1003 sage: J = random_eja()
1004 sage: x = J.random_element()
1005 sage: y = J.random_element()
1006 sage: z = J.random_element()
1007 sage: Lx = x.operator()
1008 sage: Ly = y.operator()
1009 sage: Lz = z.operator()
1010 sage: Lzy = (z*y).operator()
1011 sage: Lxy = (x*y).operator()
1012 sage: Lxz = (x*z).operator()
1013 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
1016 Test the third polarization identity from my notes or from
1019 sage: set_random_seed()
1020 sage: J = random_eja()
1021 sage: u = J.random_element()
1022 sage: y = J.random_element()
1023 sage: z = J.random_element()
1024 sage: Lu = u.operator()
1025 sage: Ly = y.operator()
1026 sage: Lz = z.operator()
1027 sage: Lzy = (z*y).operator()
1028 sage: Luy = (u*y).operator()
1029 sage: Luz = (u*z).operator()
1030 sage: Luyz = (u*(y*z)).operator()
1031 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
1032 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
1033 sage: bool(lhs == rhs)
1037 if not other
in self
.parent():
1038 raise TypeError("'other' must live in the same algebra")
1041 B
= other
.operator()
1047 Return my determinant, the product of my eigenvalues.
1051 sage: J = JordanSpinEJA(2)
1052 sage: e0,e1 = J.gens()
1053 sage: x = sum( J.gens() )
1059 sage: J = JordanSpinEJA(3)
1060 sage: e0,e1,e2 = J.gens()
1061 sage: x = sum( J.gens() )
1067 An element is invertible if and only if its determinant is
1070 sage: set_random_seed()
1071 sage: x = random_eja().random_element()
1072 sage: x.is_invertible() == (x.det() != 0)
1078 p
= P
._charpoly
_coeff
(0)
1079 # The _charpoly_coeff function already adds the factor of
1080 # -1 to ensure that _charpoly_coeff(0) is really what
1081 # appears in front of t^{0} in the charpoly. However,
1082 # we want (-1)^r times THAT for the determinant.
1083 return ((-1)**r
)*p(*self
.vector())
1088 Return the Jordan-multiplicative inverse of this element.
1092 We appeal to the quadratic representation as in Koecher's
1093 Theorem 12 in Chapter III, Section 5.
1097 The inverse in the spin factor algebra is given in Alizadeh's
1100 sage: set_random_seed()
1101 sage: n = ZZ.random_element(1,10)
1102 sage: J = JordanSpinEJA(n)
1103 sage: x = J.random_element()
1104 sage: while not x.is_invertible():
1105 ....: x = J.random_element()
1106 sage: x_vec = x.vector()
1108 sage: x_bar = x_vec[1:]
1109 sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
1110 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
1111 sage: x_inverse = coeff*inv_vec
1112 sage: x.inverse() == J(x_inverse)
1117 The identity element is its own inverse::
1119 sage: set_random_seed()
1120 sage: J = random_eja()
1121 sage: J.one().inverse() == J.one()
1124 If an element has an inverse, it acts like one::
1126 sage: set_random_seed()
1127 sage: J = random_eja()
1128 sage: x = J.random_element()
1129 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
1132 The inverse of the inverse is what we started with::
1134 sage: set_random_seed()
1135 sage: J = random_eja()
1136 sage: x = J.random_element()
1137 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
1140 The zero element is never invertible::
1142 sage: set_random_seed()
1143 sage: J = random_eja().zero().inverse()
1144 Traceback (most recent call last):
1146 ValueError: element is not invertible
1149 if not self
.is_invertible():
1150 raise ValueError("element is not invertible")
1152 return (~self
.quadratic_representation())(self
)
1155 def is_invertible(self
):
1157 Return whether or not this element is invertible.
1159 We can't use the superclass method because it relies on
1160 the algebra being associative.
1164 The usual way to do this is to check if the determinant is
1165 zero, but we need the characteristic polynomial for the
1166 determinant. The minimal polynomial is a lot easier to get,
1167 so we use Corollary 2 in Chapter V of Koecher to check
1168 whether or not the paren't algebra's zero element is a root
1169 of this element's minimal polynomial.
1173 The identity element is always invertible::
1175 sage: set_random_seed()
1176 sage: J = random_eja()
1177 sage: J.one().is_invertible()
1180 The zero element is never invertible::
1182 sage: set_random_seed()
1183 sage: J = random_eja()
1184 sage: J.zero().is_invertible()
1188 zero
= self
.parent().zero()
1189 p
= self
.minimal_polynomial()
1190 return not (p(zero
) == zero
)
1193 def is_nilpotent(self
):
1195 Return whether or not some power of this element is zero.
1197 The superclass method won't work unless we're in an
1198 associative algebra, and we aren't. However, we generate
1199 an assocoative subalgebra and we're nilpotent there if and
1200 only if we're nilpotent here (probably).
1204 The identity element is never nilpotent::
1206 sage: set_random_seed()
1207 sage: random_eja().one().is_nilpotent()
1210 The additive identity is always nilpotent::
1212 sage: set_random_seed()
1213 sage: random_eja().zero().is_nilpotent()
1217 # The element we're going to call "is_nilpotent()" on.
1218 # Either myself, interpreted as an element of a finite-
1219 # dimensional algebra, or an element of an associative
1223 if self
.parent().is_associative():
1224 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1226 V
= self
.span_of_powers()
1227 assoc_subalg
= self
.subalgebra_generated_by()
1228 # Mis-design warning: the basis used for span_of_powers()
1229 # and subalgebra_generated_by() must be the same, and in
1231 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1233 # Recursive call, but should work since elt lives in an
1234 # associative algebra.
1235 return elt
.is_nilpotent()
1238 def is_regular(self
):
1240 Return whether or not this is a regular element.
1244 The identity element always has degree one, but any element
1245 linearly-independent from it is regular::
1247 sage: J = JordanSpinEJA(5)
1248 sage: J.one().is_regular()
1250 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
1251 sage: for x in J.gens():
1252 ....: (J.one() + x).is_regular()
1260 return self
.degree() == self
.parent().rank()
1265 Compute the degree of this element the straightforward way
1266 according to the definition; by appending powers to a list
1267 and figuring out its dimension (that is, whether or not
1268 they're linearly dependent).
1272 sage: J = JordanSpinEJA(4)
1273 sage: J.one().degree()
1275 sage: e0,e1,e2,e3 = J.gens()
1276 sage: (e0 - e1).degree()
1279 In the spin factor algebra (of rank two), all elements that
1280 aren't multiples of the identity are regular::
1282 sage: set_random_seed()
1283 sage: n = ZZ.random_element(1,10)
1284 sage: J = JordanSpinEJA(n)
1285 sage: x = J.random_element()
1286 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
1290 return self
.span_of_powers().dimension()
1293 def left_matrix(self
):
1295 Our parent class defines ``left_matrix`` and ``matrix``
1296 methods whose names are misleading. We don't want them.
1298 raise NotImplementedError("use operator().matrix() instead")
1300 matrix
= left_matrix
1303 def minimal_polynomial(self
):
1305 Return the minimal polynomial of this element,
1306 as a function of the variable `t`.
1310 We restrict ourselves to the associative subalgebra
1311 generated by this element, and then return the minimal
1312 polynomial of this element's operator matrix (in that
1313 subalgebra). This works by Baes Proposition 2.3.16.
1317 The minimal polynomial of the identity and zero elements are
1320 sage: set_random_seed()
1321 sage: J = random_eja()
1322 sage: J.one().minimal_polynomial()
1324 sage: J.zero().minimal_polynomial()
1327 The degree of an element is (by one definition) the degree
1328 of its minimal polynomial::
1330 sage: set_random_seed()
1331 sage: x = random_eja().random_element()
1332 sage: x.degree() == x.minimal_polynomial().degree()
1335 The minimal polynomial and the characteristic polynomial coincide
1336 and are known (see Alizadeh, Example 11.11) for all elements of
1337 the spin factor algebra that aren't scalar multiples of the
1340 sage: set_random_seed()
1341 sage: n = ZZ.random_element(2,10)
1342 sage: J = JordanSpinEJA(n)
1343 sage: y = J.random_element()
1344 sage: while y == y.coefficient(0)*J.one():
1345 ....: y = J.random_element()
1346 sage: y0 = y.vector()[0]
1347 sage: y_bar = y.vector()[1:]
1348 sage: actual = y.minimal_polynomial()
1349 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1350 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1351 sage: bool(actual == expected)
1354 The minimal polynomial should always kill its element::
1356 sage: set_random_seed()
1357 sage: x = random_eja().random_element()
1358 sage: p = x.minimal_polynomial()
1359 sage: x.apply_univariate_polynomial(p)
1363 V
= self
.span_of_powers()
1364 assoc_subalg
= self
.subalgebra_generated_by()
1365 # Mis-design warning: the basis used for span_of_powers()
1366 # and subalgebra_generated_by() must be the same, and in
1368 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1369 return elt
.operator().minimal_polynomial()
1373 def natural_representation(self
):
1375 Return a more-natural representation of this element.
1377 Every finite-dimensional Euclidean Jordan Algebra is a
1378 direct sum of five simple algebras, four of which comprise
1379 Hermitian matrices. This method returns the original
1380 "natural" representation of this element as a Hermitian
1381 matrix, if it has one. If not, you get the usual representation.
1385 sage: J = ComplexHermitianEJA(3)
1388 sage: J.one().natural_representation()
1398 sage: J = QuaternionHermitianEJA(3)
1401 sage: J.one().natural_representation()
1402 [1 0 0 0 0 0 0 0 0 0 0 0]
1403 [0 1 0 0 0 0 0 0 0 0 0 0]
1404 [0 0 1 0 0 0 0 0 0 0 0 0]
1405 [0 0 0 1 0 0 0 0 0 0 0 0]
1406 [0 0 0 0 1 0 0 0 0 0 0 0]
1407 [0 0 0 0 0 1 0 0 0 0 0 0]
1408 [0 0 0 0 0 0 1 0 0 0 0 0]
1409 [0 0 0 0 0 0 0 1 0 0 0 0]
1410 [0 0 0 0 0 0 0 0 1 0 0 0]
1411 [0 0 0 0 0 0 0 0 0 1 0 0]
1412 [0 0 0 0 0 0 0 0 0 0 1 0]
1413 [0 0 0 0 0 0 0 0 0 0 0 1]
1416 B
= self
.parent().natural_basis()
1417 W
= B
[0].matrix_space()
1418 return W
.linear_combination(zip(self
.vector(), B
))
1423 Return the left-multiplication-by-this-element
1424 operator on the ambient algebra.
1428 sage: set_random_seed()
1429 sage: J = random_eja()
1430 sage: x = J.random_element()
1431 sage: y = J.random_element()
1432 sage: x.operator()(y) == x*y
1434 sage: y.operator()(x) == x*y
1439 fda_elt
= FiniteDimensionalAlgebraElement(P
, self
)
1440 return FiniteDimensionalEuclideanJordanAlgebraOperator(
1443 fda_elt
.matrix().transpose() )
1446 def quadratic_representation(self
, other
=None):
1448 Return the quadratic representation of this element.
1452 The explicit form in the spin factor algebra is given by
1453 Alizadeh's Example 11.12::
1455 sage: set_random_seed()
1456 sage: n = ZZ.random_element(1,10)
1457 sage: J = JordanSpinEJA(n)
1458 sage: x = J.random_element()
1459 sage: x_vec = x.vector()
1461 sage: x_bar = x_vec[1:]
1462 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1463 sage: B = 2*x0*x_bar.row()
1464 sage: C = 2*x0*x_bar.column()
1465 sage: D = identity_matrix(QQ, n-1)
1466 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1467 sage: D = D + 2*x_bar.tensor_product(x_bar)
1468 sage: Q = block_matrix(2,2,[A,B,C,D])
1469 sage: Q == x.quadratic_representation().matrix()
1472 Test all of the properties from Theorem 11.2 in Alizadeh::
1474 sage: set_random_seed()
1475 sage: J = random_eja()
1476 sage: x = J.random_element()
1477 sage: y = J.random_element()
1478 sage: Lx = x.operator()
1479 sage: Lxx = (x*x).operator()
1480 sage: Qx = x.quadratic_representation()
1481 sage: Qy = y.quadratic_representation()
1482 sage: Qxy = x.quadratic_representation(y)
1483 sage: Qex = J.one().quadratic_representation(x)
1484 sage: n = ZZ.random_element(10)
1485 sage: Qxn = (x^n).quadratic_representation()
1489 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1492 Property 2 (multiply on the right for :trac:`28272`):
1494 sage: alpha = QQ.random_element()
1495 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1500 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1503 sage: not x.is_invertible() or (
1506 ....: x.inverse().quadratic_representation() )
1509 sage: Qxy(J.one()) == x*y
1514 sage: not x.is_invertible() or (
1515 ....: x.quadratic_representation(x.inverse())*Qx
1516 ....: == Qx*x.quadratic_representation(x.inverse()) )
1519 sage: not x.is_invertible() or (
1520 ....: x.quadratic_representation(x.inverse())*Qx
1522 ....: 2*x.operator()*Qex - Qx )
1525 sage: 2*x.operator()*Qex - Qx == Lxx
1530 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1540 sage: not x.is_invertible() or (
1541 ....: Qx*x.inverse().operator() == Lx )
1546 sage: not x.operator_commutes_with(y) or (
1547 ....: Qx(y)^n == Qxn(y^n) )
1553 elif not other
in self
.parent():
1554 raise TypeError("'other' must live in the same algebra")
1557 M
= other
.operator()
1558 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1561 def span_of_powers(self
):
1563 Return the vector space spanned by successive powers of
1566 # The dimension of the subalgebra can't be greater than
1567 # the big algebra, so just put everything into a list
1568 # and let span() get rid of the excess.
1570 # We do the extra ambient_vector_space() in case we're messing
1571 # with polynomials and the direct parent is a module.
1572 V
= self
.parent().vector_space()
1573 return V
.span( (self
**d
).vector() for d
in xrange(V
.dimension()) )
1576 def subalgebra_generated_by(self
):
1578 Return the associative subalgebra of the parent EJA generated
1583 sage: set_random_seed()
1584 sage: x = random_eja().random_element()
1585 sage: x.subalgebra_generated_by().is_associative()
1588 Squaring in the subalgebra should work the same as in
1591 sage: set_random_seed()
1592 sage: x = random_eja().random_element()
1593 sage: u = x.subalgebra_generated_by().random_element()
1594 sage: u.operator()(u) == u^2
1598 # First get the subspace spanned by the powers of myself...
1599 V
= self
.span_of_powers()
1600 F
= self
.base_ring()
1602 # Now figure out the entries of the right-multiplication
1603 # matrix for the successive basis elements b0, b1,... of
1606 for b_right
in V
.basis():
1607 eja_b_right
= self
.parent()(b_right
)
1609 # The first row of the right-multiplication matrix by
1610 # b1 is what we get if we apply that matrix to b1. The
1611 # second row of the right multiplication matrix by b1
1612 # is what we get when we apply that matrix to b2...
1614 # IMPORTANT: this assumes that all vectors are COLUMN
1615 # vectors, unlike our superclass (which uses row vectors).
1616 for b_left
in V
.basis():
1617 eja_b_left
= self
.parent()(b_left
)
1618 # Multiply in the original EJA, but then get the
1619 # coordinates from the subalgebra in terms of its
1621 this_row
= V
.coordinates((eja_b_left
*eja_b_right
).vector())
1622 b_right_rows
.append(this_row
)
1623 b_right_matrix
= matrix(F
, b_right_rows
)
1624 mats
.append(b_right_matrix
)
1626 # It's an algebra of polynomials in one element, and EJAs
1627 # are power-associative.
1629 # TODO: choose generator names intelligently.
1630 return FiniteDimensionalEuclideanJordanAlgebra(F
, mats
, assume_associative
=True, names
='f')
1633 def subalgebra_idempotent(self
):
1635 Find an idempotent in the associative subalgebra I generate
1636 using Proposition 2.3.5 in Baes.
1640 sage: set_random_seed()
1641 sage: J = random_eja()
1642 sage: x = J.random_element()
1643 sage: while x.is_nilpotent():
1644 ....: x = J.random_element()
1645 sage: c = x.subalgebra_idempotent()
1650 if self
.is_nilpotent():
1651 raise ValueError("this only works with non-nilpotent elements!")
1653 V
= self
.span_of_powers()
1654 J
= self
.subalgebra_generated_by()
1655 # Mis-design warning: the basis used for span_of_powers()
1656 # and subalgebra_generated_by() must be the same, and in
1658 u
= J(V
.coordinates(self
.vector()))
1660 # The image of the matrix of left-u^m-multiplication
1661 # will be minimal for some natural number s...
1663 minimal_dim
= V
.dimension()
1664 for i
in xrange(1, V
.dimension()):
1665 this_dim
= (u
**i
).operator().matrix().image().dimension()
1666 if this_dim
< minimal_dim
:
1667 minimal_dim
= this_dim
1670 # Now minimal_matrix should correspond to the smallest
1671 # non-zero subspace in Baes's (or really, Koecher's)
1674 # However, we need to restrict the matrix to work on the
1675 # subspace... or do we? Can't we just solve, knowing that
1676 # A(c) = u^(s+1) should have a solution in the big space,
1679 # Beware, solve_right() means that we're using COLUMN vectors.
1680 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1682 A
= u_next
.operator().matrix()
1683 c_coordinates
= A
.solve_right(u_next
.vector())
1685 # Now c_coordinates is the idempotent we want, but it's in
1686 # the coordinate system of the subalgebra.
1688 # We need the basis for J, but as elements of the parent algebra.
1690 basis
= [self
.parent(v
) for v
in V
.basis()]
1691 return self
.parent().linear_combination(zip(c_coordinates
, basis
))
1696 Return my trace, the sum of my eigenvalues.
1700 sage: J = JordanSpinEJA(3)
1701 sage: x = sum(J.gens())
1707 sage: J = RealCartesianProductEJA(5)
1708 sage: J.one().trace()
1713 The trace of an element is a real number::
1715 sage: set_random_seed()
1716 sage: J = random_eja()
1717 sage: J.random_element().trace() in J.base_ring()
1723 p
= P
._charpoly
_coeff
(r
-1)
1724 # The _charpoly_coeff function already adds the factor of
1725 # -1 to ensure that _charpoly_coeff(r-1) is really what
1726 # appears in front of t^{r-1} in the charpoly. However,
1727 # we want the negative of THAT for the trace.
1728 return -p(*self
.vector())
1731 def trace_inner_product(self
, other
):
1733 Return the trace inner product of myself and ``other``.
1737 The trace inner product is commutative::
1739 sage: set_random_seed()
1740 sage: J = random_eja()
1741 sage: x = J.random_element(); y = J.random_element()
1742 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1745 The trace inner product is bilinear::
1747 sage: set_random_seed()
1748 sage: J = random_eja()
1749 sage: x = J.random_element()
1750 sage: y = J.random_element()
1751 sage: z = J.random_element()
1752 sage: a = QQ.random_element();
1753 sage: actual = (a*(x+z)).trace_inner_product(y)
1754 sage: expected = ( a*x.trace_inner_product(y) +
1755 ....: a*z.trace_inner_product(y) )
1756 sage: actual == expected
1758 sage: actual = x.trace_inner_product(a*(y+z))
1759 sage: expected = ( a*x.trace_inner_product(y) +
1760 ....: a*x.trace_inner_product(z) )
1761 sage: actual == expected
1764 The trace inner product satisfies the compatibility
1765 condition in the definition of a Euclidean Jordan algebra::
1767 sage: set_random_seed()
1768 sage: J = random_eja()
1769 sage: x = J.random_element()
1770 sage: y = J.random_element()
1771 sage: z = J.random_element()
1772 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1776 if not other
in self
.parent():
1777 raise TypeError("'other' must live in the same algebra")
1779 return (self
*other
).trace()
1782 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1784 Return the Euclidean Jordan Algebra corresponding to the set
1785 `R^n` under the Hadamard product.
1787 Note: this is nothing more than the Cartesian product of ``n``
1788 copies of the spin algebra. Once Cartesian product algebras
1789 are implemented, this can go.
1793 This multiplication table can be verified by hand::
1795 sage: J = RealCartesianProductEJA(3)
1796 sage: e0,e1,e2 = J.gens()
1812 def __classcall_private__(cls
, n
, field
=QQ
):
1813 # The FiniteDimensionalAlgebra constructor takes a list of
1814 # matrices, the ith representing right multiplication by the ith
1815 # basis element in the vector space. So if e_1 = (1,0,0), then
1816 # right (Hadamard) multiplication of x by e_1 picks out the first
1817 # component of x; and likewise for the ith basis element e_i.
1818 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
1819 for i
in xrange(n
) ]
1821 fdeja
= super(RealCartesianProductEJA
, cls
)
1822 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=n
)
1824 def inner_product(self
, x
, y
):
1825 return _usual_ip(x
,y
)
1830 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1834 For now, we choose a random natural number ``n`` (greater than zero)
1835 and then give you back one of the following:
1837 * The cartesian product of the rational numbers ``n`` times; this is
1838 ``QQ^n`` with the Hadamard product.
1840 * The Jordan spin algebra on ``QQ^n``.
1842 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1845 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1846 in the space of ``2n``-by-``2n`` real symmetric matrices.
1848 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1849 in the space of ``4n``-by-``4n`` real symmetric matrices.
1851 Later this might be extended to return Cartesian products of the
1857 Euclidean Jordan algebra of degree...
1861 # The max_n component lets us choose different upper bounds on the
1862 # value "n" that gets passed to the constructor. This is needed
1863 # because e.g. R^{10} is reasonable to test, while the Hermitian
1864 # 10-by-10 quaternion matrices are not.
1865 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
1867 (RealSymmetricEJA
, 5),
1868 (ComplexHermitianEJA
, 4),
1869 (QuaternionHermitianEJA
, 3)])
1870 n
= ZZ
.random_element(1, max_n
)
1871 return constructor(n
, field
=QQ
)
1875 def _real_symmetric_basis(n
, field
=QQ
):
1877 Return a basis for the space of real symmetric n-by-n matrices.
1879 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1883 for j
in xrange(i
+1):
1884 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1888 # Beware, orthogonal but not normalized!
1889 Sij
= Eij
+ Eij
.transpose()
1894 def _complex_hermitian_basis(n
, field
=QQ
):
1896 Returns a basis for the space of complex Hermitian n-by-n matrices.
1900 sage: set_random_seed()
1901 sage: n = ZZ.random_element(1,5)
1902 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1906 F
= QuadraticField(-1, 'I')
1909 # This is like the symmetric case, but we need to be careful:
1911 # * We want conjugate-symmetry, not just symmetry.
1912 # * The diagonal will (as a result) be real.
1916 for j
in xrange(i
+1):
1917 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1919 Sij
= _embed_complex_matrix(Eij
)
1922 # Beware, orthogonal but not normalized! The second one
1923 # has a minus because it's conjugated.
1924 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
1926 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
1931 def _quaternion_hermitian_basis(n
, field
=QQ
):
1933 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1937 sage: set_random_seed()
1938 sage: n = ZZ.random_element(1,5)
1939 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1943 Q
= QuaternionAlgebra(QQ
,-1,-1)
1946 # This is like the symmetric case, but we need to be careful:
1948 # * We want conjugate-symmetry, not just symmetry.
1949 # * The diagonal will (as a result) be real.
1953 for j
in xrange(i
+1):
1954 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1956 Sij
= _embed_quaternion_matrix(Eij
)
1959 # Beware, orthogonal but not normalized! The second,
1960 # third, and fourth ones have a minus because they're
1962 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
1964 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
1966 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
1968 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
1974 return vector(m
.base_ring(), m
.list())
1977 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
1979 def _multiplication_table_from_matrix_basis(basis
):
1981 At least three of the five simple Euclidean Jordan algebras have the
1982 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1983 multiplication on the right is matrix multiplication. Given a basis
1984 for the underlying matrix space, this function returns a
1985 multiplication table (obtained by looping through the basis
1986 elements) for an algebra of those matrices. A reordered copy
1987 of the basis is also returned to work around the fact that
1988 the ``span()`` in this function will change the order of the basis
1989 from what we think it is, to... something else.
1991 # In S^2, for example, we nominally have four coordinates even
1992 # though the space is of dimension three only. The vector space V
1993 # is supposed to hold the entire long vector, and the subspace W
1994 # of V will be spanned by the vectors that arise from symmetric
1995 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1996 field
= basis
[0].base_ring()
1997 dimension
= basis
[0].nrows()
1999 V
= VectorSpace(field
, dimension
**2)
2000 W
= V
.span( _mat2vec(s
) for s
in basis
)
2002 # Taking the span above reorders our basis (thanks, jerk!) so we
2003 # need to put our "matrix basis" in the same order as the
2004 # (reordered) vector basis.
2005 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
2009 # Brute force the multiplication-by-s matrix by looping
2010 # through all elements of the basis and doing the computation
2011 # to find out what the corresponding row should be. BEWARE:
2012 # these multiplication tables won't be symmetric! It therefore
2013 # becomes REALLY IMPORTANT that the underlying algebra
2014 # constructor uses ROW vectors and not COLUMN vectors. That's
2015 # why we're computing rows here and not columns.
2018 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
2019 Q_rows
.append(W
.coordinates(this_row
))
2020 Q
= matrix(field
, W
.dimension(), Q_rows
)
2026 def _embed_complex_matrix(M
):
2028 Embed the n-by-n complex matrix ``M`` into the space of real
2029 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2030 bi` to the block matrix ``[[a,b],[-b,a]]``.
2034 sage: F = QuadraticField(-1,'i')
2035 sage: x1 = F(4 - 2*i)
2036 sage: x2 = F(1 + 2*i)
2039 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2040 sage: _embed_complex_matrix(M)
2049 Embedding is a homomorphism (isomorphism, in fact)::
2051 sage: set_random_seed()
2052 sage: n = ZZ.random_element(5)
2053 sage: F = QuadraticField(-1, 'i')
2054 sage: X = random_matrix(F, n)
2055 sage: Y = random_matrix(F, n)
2056 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
2057 sage: expected = _embed_complex_matrix(X*Y)
2058 sage: actual == expected
2064 raise ValueError("the matrix 'M' must be square")
2065 field
= M
.base_ring()
2070 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
2072 # We can drop the imaginaries here.
2073 return block_matrix(field
.base_ring(), n
, blocks
)
2076 def _unembed_complex_matrix(M
):
2078 The inverse of _embed_complex_matrix().
2082 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2083 ....: [-2, 1, -4, 3],
2084 ....: [ 9, 10, 11, 12],
2085 ....: [-10, 9, -12, 11] ])
2086 sage: _unembed_complex_matrix(A)
2088 [ 10*i + 9 12*i + 11]
2092 Unembedding is the inverse of embedding::
2094 sage: set_random_seed()
2095 sage: F = QuadraticField(-1, 'i')
2096 sage: M = random_matrix(F, 3)
2097 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
2103 raise ValueError("the matrix 'M' must be square")
2104 if not n
.mod(2).is_zero():
2105 raise ValueError("the matrix 'M' must be a complex embedding")
2107 F
= QuadraticField(-1, 'i')
2110 # Go top-left to bottom-right (reading order), converting every
2111 # 2-by-2 block we see to a single complex element.
2113 for k
in xrange(n
/2):
2114 for j
in xrange(n
/2):
2115 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
2116 if submat
[0,0] != submat
[1,1]:
2117 raise ValueError('bad on-diagonal submatrix')
2118 if submat
[0,1] != -submat
[1,0]:
2119 raise ValueError('bad off-diagonal submatrix')
2120 z
= submat
[0,0] + submat
[0,1]*i
2123 return matrix(F
, n
/2, elements
)
2126 def _embed_quaternion_matrix(M
):
2128 Embed the n-by-n quaternion matrix ``M`` into the space of real
2129 matrices of size 4n-by-4n by first sending each quaternion entry
2130 `z = a + bi + cj + dk` to the block-complex matrix
2131 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
2136 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2137 sage: i,j,k = Q.gens()
2138 sage: x = 1 + 2*i + 3*j + 4*k
2139 sage: M = matrix(Q, 1, [[x]])
2140 sage: _embed_quaternion_matrix(M)
2146 Embedding is a homomorphism (isomorphism, in fact)::
2148 sage: set_random_seed()
2149 sage: n = ZZ.random_element(5)
2150 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2151 sage: X = random_matrix(Q, n)
2152 sage: Y = random_matrix(Q, n)
2153 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
2154 sage: expected = _embed_quaternion_matrix(X*Y)
2155 sage: actual == expected
2159 quaternions
= M
.base_ring()
2162 raise ValueError("the matrix 'M' must be square")
2164 F
= QuadraticField(-1, 'i')
2169 t
= z
.coefficient_tuple()
2174 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2175 [-c
+ d
*i
, a
- b
*i
]])
2176 blocks
.append(_embed_complex_matrix(cplx_matrix
))
2178 # We should have real entries by now, so use the realest field
2179 # we've got for the return value.
2180 return block_matrix(quaternions
.base_ring(), n
, blocks
)
2183 def _unembed_quaternion_matrix(M
):
2185 The inverse of _embed_quaternion_matrix().
2189 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2190 ....: [-2, 1, -4, 3],
2191 ....: [-3, 4, 1, -2],
2192 ....: [-4, -3, 2, 1]])
2193 sage: _unembed_quaternion_matrix(M)
2194 [1 + 2*i + 3*j + 4*k]
2198 Unembedding is the inverse of embedding::
2200 sage: set_random_seed()
2201 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2202 sage: M = random_matrix(Q, 3)
2203 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
2209 raise ValueError("the matrix 'M' must be square")
2210 if not n
.mod(4).is_zero():
2211 raise ValueError("the matrix 'M' must be a complex embedding")
2213 Q
= QuaternionAlgebra(QQ
,-1,-1)
2216 # Go top-left to bottom-right (reading order), converting every
2217 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2220 for l
in xrange(n
/4):
2221 for m
in xrange(n
/4):
2222 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
2223 if submat
[0,0] != submat
[1,1].conjugate():
2224 raise ValueError('bad on-diagonal submatrix')
2225 if submat
[0,1] != -submat
[1,0].conjugate():
2226 raise ValueError('bad off-diagonal submatrix')
2227 z
= submat
[0,0].real() + submat
[0,0].imag()*i
2228 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
2231 return matrix(Q
, n
/4, elements
)
2234 # The usual inner product on R^n.
2236 return x
.vector().inner_product(y
.vector())
2238 # The inner product used for the real symmetric simple EJA.
2239 # We keep it as a separate function because e.g. the complex
2240 # algebra uses the same inner product, except divided by 2.
2241 def _matrix_ip(X
,Y
):
2242 X_mat
= X
.natural_representation()
2243 Y_mat
= Y
.natural_representation()
2244 return (X_mat
*Y_mat
).trace()
2247 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2249 The rank-n simple EJA consisting of real symmetric n-by-n
2250 matrices, the usual symmetric Jordan product, and the trace inner
2251 product. It has dimension `(n^2 + n)/2` over the reals.
2255 sage: J = RealSymmetricEJA(2)
2256 sage: e0, e1, e2 = J.gens()
2266 The degree of this algebra is `(n^2 + n) / 2`::
2268 sage: set_random_seed()
2269 sage: n = ZZ.random_element(1,5)
2270 sage: J = RealSymmetricEJA(n)
2271 sage: J.degree() == (n^2 + n)/2
2274 The Jordan multiplication is what we think it is::
2276 sage: set_random_seed()
2277 sage: n = ZZ.random_element(1,5)
2278 sage: J = RealSymmetricEJA(n)
2279 sage: x = J.random_element()
2280 sage: y = J.random_element()
2281 sage: actual = (x*y).natural_representation()
2282 sage: X = x.natural_representation()
2283 sage: Y = y.natural_representation()
2284 sage: expected = (X*Y + Y*X)/2
2285 sage: actual == expected
2287 sage: J(expected) == x*y
2292 def __classcall_private__(cls
, n
, field
=QQ
):
2293 S
= _real_symmetric_basis(n
, field
=field
)
2294 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2296 fdeja
= super(RealSymmetricEJA
, cls
)
2297 return fdeja
.__classcall
_private
__(cls
,
2303 def inner_product(self
, x
, y
):
2304 return _matrix_ip(x
,y
)
2307 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2309 The rank-n simple EJA consisting of complex Hermitian n-by-n
2310 matrices over the real numbers, the usual symmetric Jordan product,
2311 and the real-part-of-trace inner product. It has dimension `n^2` over
2316 The degree of this algebra is `n^2`::
2318 sage: set_random_seed()
2319 sage: n = ZZ.random_element(1,5)
2320 sage: J = ComplexHermitianEJA(n)
2321 sage: J.degree() == n^2
2324 The Jordan multiplication is what we think it is::
2326 sage: set_random_seed()
2327 sage: n = ZZ.random_element(1,5)
2328 sage: J = ComplexHermitianEJA(n)
2329 sage: x = J.random_element()
2330 sage: y = J.random_element()
2331 sage: actual = (x*y).natural_representation()
2332 sage: X = x.natural_representation()
2333 sage: Y = y.natural_representation()
2334 sage: expected = (X*Y + Y*X)/2
2335 sage: actual == expected
2337 sage: J(expected) == x*y
2342 def __classcall_private__(cls
, n
, field
=QQ
):
2343 S
= _complex_hermitian_basis(n
)
2344 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2346 fdeja
= super(ComplexHermitianEJA
, cls
)
2347 return fdeja
.__classcall
_private
__(cls
,
2353 def inner_product(self
, x
, y
):
2354 # Since a+bi on the diagonal is represented as
2359 # we'll double-count the "a" entries if we take the trace of
2361 return _matrix_ip(x
,y
)/2
2364 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2366 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2367 matrices, the usual symmetric Jordan product, and the
2368 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2373 The degree of this algebra is `n^2`::
2375 sage: set_random_seed()
2376 sage: n = ZZ.random_element(1,5)
2377 sage: J = QuaternionHermitianEJA(n)
2378 sage: J.degree() == 2*(n^2) - n
2381 The Jordan multiplication is what we think it is::
2383 sage: set_random_seed()
2384 sage: n = ZZ.random_element(1,5)
2385 sage: J = QuaternionHermitianEJA(n)
2386 sage: x = J.random_element()
2387 sage: y = J.random_element()
2388 sage: actual = (x*y).natural_representation()
2389 sage: X = x.natural_representation()
2390 sage: Y = y.natural_representation()
2391 sage: expected = (X*Y + Y*X)/2
2392 sage: actual == expected
2394 sage: J(expected) == x*y
2399 def __classcall_private__(cls
, n
, field
=QQ
):
2400 S
= _quaternion_hermitian_basis(n
)
2401 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2403 fdeja
= super(QuaternionHermitianEJA
, cls
)
2404 return fdeja
.__classcall
_private
__(cls
,
2410 def inner_product(self
, x
, y
):
2411 # Since a+bi+cj+dk on the diagonal is represented as
2413 # a + bi +cj + dk = [ a b c d]
2418 # we'll quadruple-count the "a" entries if we take the trace of
2420 return _matrix_ip(x
,y
)/4
2423 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2425 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2426 with the usual inner product and jordan product ``x*y =
2427 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2432 This multiplication table can be verified by hand::
2434 sage: J = JordanSpinEJA(4)
2435 sage: e0,e1,e2,e3 = J.gens()
2453 def __classcall_private__(cls
, n
, field
=QQ
):
2455 id_matrix
= identity_matrix(field
, n
)
2457 ei
= id_matrix
.column(i
)
2458 Qi
= zero_matrix(field
, n
)
2460 Qi
.set_column(0, ei
)
2461 Qi
+= diagonal_matrix(n
, [ei
[0]]*n
)
2462 # The addition of the diagonal matrix adds an extra ei[0] in the
2463 # upper-left corner of the matrix.
2464 Qi
[0,0] = Qi
[0,0] * ~
field(2)
2467 # The rank of the spin algebra is two, unless we're in a
2468 # one-dimensional ambient space (because the rank is bounded by
2469 # the ambient dimension).
2470 fdeja
= super(JordanSpinEJA
, cls
)
2471 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=min(n
,2))
2473 def inner_product(self
, x
, y
):
2474 return _usual_ip(x
,y
)