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 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra
):
356 def __classcall_private__(cls
,
360 assume_associative
=False,
365 mult_table
= [b
.base_extend(field
) for b
in mult_table
]
368 if not (is_Matrix(b
) and b
.dimensions() == (n
, n
)):
369 raise ValueError("input is not a multiplication table")
370 mult_table
= tuple(mult_table
)
372 cat
= FiniteDimensionalAlgebrasWithBasis(field
)
373 cat
.or_subcategory(category
)
374 if assume_associative
:
375 cat
= cat
.Associative()
377 names
= normalize_names(n
, names
)
379 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, cls
)
380 return fda
.__classcall
__(cls
,
383 assume_associative
=assume_associative
,
387 natural_basis
=natural_basis
)
394 assume_associative
=False,
401 By definition, Jordan multiplication commutes::
403 sage: set_random_seed()
404 sage: J = random_eja()
405 sage: x = J.random_element()
406 sage: y = J.random_element()
412 self
._natural
_basis
= natural_basis
413 self
._multiplication
_table
= mult_table
414 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
423 Return a string representation of ``self``.
425 fmt
= "Euclidean Jordan algebra of degree {} over {}"
426 return fmt
.format(self
.degree(), self
.base_ring())
429 def _a_regular_element(self
):
431 Guess a regular element. Needed to compute the basis for our
432 characteristic polynomial coefficients.
435 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
436 if not z
.is_regular():
437 raise ValueError("don't know a regular element")
442 def _charpoly_basis_space(self
):
444 Return the vector space spanned by the basis used in our
445 characteristic polynomial coefficients. This is used not only to
446 compute those coefficients, but also any time we need to
447 evaluate the coefficients (like when we compute the trace or
450 z
= self
._a
_regular
_element
()
451 V
= self
.vector_space()
452 V1
= V
.span_of_basis( (z
**k
).vector() for k
in range(self
.rank()) )
453 b
= (V1
.basis() + V1
.complement().basis())
454 return V
.span_of_basis(b
)
458 def _charpoly_coeff(self
, i
):
460 Return the coefficient polynomial "a_{i}" of this algebra's
461 general characteristic polynomial.
463 Having this be a separate cached method lets us compute and
464 store the trace/determinant (a_{r-1} and a_{0} respectively)
465 separate from the entire characteristic polynomial.
467 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
468 R
= A_of_x
.base_ring()
470 # Guaranteed by theory
473 # Danger: the in-place modification is done for performance
474 # reasons (reconstructing a matrix with huge polynomial
475 # entries is slow), but I don't know how cached_method works,
476 # so it's highly possible that we're modifying some global
477 # list variable by reference, here. In other words, you
478 # probably shouldn't call this method twice on the same
479 # algebra, at the same time, in two threads
480 Ai_orig
= A_of_x
.column(i
)
481 A_of_x
.set_column(i
,xr
)
482 numerator
= A_of_x
.det()
483 A_of_x
.set_column(i
,Ai_orig
)
485 # We're relying on the theory here to ensure that each a_i is
486 # indeed back in R, and the added negative signs are to make
487 # the whole charpoly expression sum to zero.
488 return R(-numerator
/detA
)
492 def _charpoly_matrix_system(self
):
494 Compute the matrix whose entries A_ij are polynomials in
495 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
496 corresponding to `x^r` and the determinent of the matrix A =
497 [A_ij]. In other words, all of the fixed (cachable) data needed
498 to compute the coefficients of the characteristic polynomial.
503 # Construct a new algebra over a multivariate polynomial ring...
504 names
= ['X' + str(i
) for i
in range(1,n
+1)]
505 R
= PolynomialRing(self
.base_ring(), names
)
506 J
= FiniteDimensionalEuclideanJordanAlgebra(R
,
507 self
._multiplication
_table
,
510 idmat
= identity_matrix(J
.base_ring(), n
)
512 W
= self
._charpoly
_basis
_space
()
513 W
= W
.change_ring(R
.fraction_field())
515 # Starting with the standard coordinates x = (X1,X2,...,Xn)
516 # and then converting the entries to W-coordinates allows us
517 # to pass in the standard coordinates to the charpoly and get
518 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
521 # W.coordinates(x^2) eval'd at (standard z-coords)
525 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
527 # We want the middle equivalent thing in our matrix, but use
528 # the first equivalent thing instead so that we can pass in
529 # standard coordinates.
530 x
= J(vector(R
, R
.gens()))
531 l1
= [column_matrix(W
.coordinates((x
**k
).vector())) for k
in range(r
)]
532 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
533 A_of_x
= block_matrix(R
, 1, n
, (l1
+ l2
))
534 xr
= W
.coordinates((x
**r
).vector())
535 return (A_of_x
, x
, xr
, A_of_x
.det())
539 def characteristic_polynomial(self
):
544 This implementation doesn't guarantee that the polynomial
545 denominator in the coefficients is not identically zero, so
546 theoretically it could crash. The way that this is handled
547 in e.g. Faraut and Koranyi is to use a basis that guarantees
548 the denominator is non-zero. But, doing so requires knowledge
549 of at least one regular element, and we don't even know how
550 to do that. The trade-off is that, if we use the standard basis,
551 the resulting polynomial will accept the "usual" coordinates. In
552 other words, we don't have to do a change of basis before e.g.
553 computing the trace or determinant.
557 The characteristic polynomial in the spin algebra is given in
558 Alizadeh, Example 11.11::
560 sage: J = JordanSpinEJA(3)
561 sage: p = J.characteristic_polynomial(); p
562 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
563 sage: xvec = J.one().vector()
571 # The list of coefficient polynomials a_1, a_2, ..., a_n.
572 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
574 # We go to a bit of trouble here to reorder the
575 # indeterminates, so that it's easier to evaluate the
576 # characteristic polynomial at x's coordinates and get back
577 # something in terms of t, which is what we want.
579 S
= PolynomialRing(self
.base_ring(),'t')
581 S
= PolynomialRing(S
, R
.variable_names())
584 # Note: all entries past the rth should be zero. The
585 # coefficient of the highest power (x^r) is 1, but it doesn't
586 # appear in the solution vector which contains coefficients
587 # for the other powers (to make them sum to x^r).
589 a
[r
] = 1 # corresponds to x^r
591 # When the rank is equal to the dimension, trying to
592 # assign a[r] goes out-of-bounds.
593 a
.append(1) # corresponds to x^r
595 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
598 def inner_product(self
, x
, y
):
600 The inner product associated with this Euclidean Jordan algebra.
602 Defaults to the trace inner product, but can be overridden by
603 subclasses if they are sure that the necessary properties are
608 The inner product must satisfy its axiom for this algebra to truly
609 be a Euclidean Jordan Algebra::
611 sage: set_random_seed()
612 sage: J = random_eja()
613 sage: x = J.random_element()
614 sage: y = J.random_element()
615 sage: z = J.random_element()
616 sage: (x*y).inner_product(z) == y.inner_product(x*z)
620 if (not x
in self
) or (not y
in self
):
621 raise TypeError("arguments must live in this algebra")
622 return x
.trace_inner_product(y
)
625 def natural_basis(self
):
627 Return a more-natural representation of this algebra's basis.
629 Every finite-dimensional Euclidean Jordan Algebra is a direct
630 sum of five simple algebras, four of which comprise Hermitian
631 matrices. This method returns the original "natural" basis
632 for our underlying vector space. (Typically, the natural basis
633 is used to construct the multiplication table in the first place.)
635 Note that this will always return a matrix. The standard basis
636 in `R^n` will be returned as `n`-by-`1` column matrices.
640 sage: J = RealSymmetricEJA(2)
643 sage: J.natural_basis()
651 sage: J = JordanSpinEJA(2)
654 sage: J.natural_basis()
661 if self
._natural
_basis
is None:
662 return tuple( b
.vector().column() for b
in self
.basis() )
664 return self
._natural
_basis
669 Return the rank of this EJA.
671 if self
._rank
is None:
672 raise ValueError("no rank specified at genesis")
676 def vector_space(self
):
678 Return the vector space that underlies this algebra.
682 sage: J = RealSymmetricEJA(2)
683 sage: J.vector_space()
684 Vector space of dimension 3 over Rational Field
687 return self
.zero().vector().parent().ambient_vector_space()
690 class Element(FiniteDimensionalAlgebraElement
):
692 An element of a Euclidean Jordan algebra.
697 Oh man, I should not be doing this. This hides the "disabled"
698 methods ``left_matrix`` and ``matrix`` from introspection;
699 in particular it removes them from tab-completion.
701 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
702 dir(self
.__class
__) )
705 def __init__(self
, A
, elt
=None):
709 The identity in `S^n` is converted to the identity in the EJA::
711 sage: J = RealSymmetricEJA(3)
712 sage: I = identity_matrix(QQ,3)
713 sage: J(I) == J.one()
716 This skew-symmetric matrix can't be represented in the EJA::
718 sage: J = RealSymmetricEJA(3)
719 sage: A = matrix(QQ,3, lambda i,j: i-j)
721 Traceback (most recent call last):
723 ArithmeticError: vector is not in free module
726 # Goal: if we're given a matrix, and if it lives in our
727 # parent algebra's "natural ambient space," convert it
728 # into an algebra element.
730 # The catch is, we make a recursive call after converting
731 # the given matrix into a vector that lives in the algebra.
732 # This we need to try the parent class initializer first,
733 # to avoid recursing forever if we're given something that
734 # already fits into the algebra, but also happens to live
735 # in the parent's "natural ambient space" (this happens with
738 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, elt
)
740 natural_basis
= A
.natural_basis()
741 if elt
in natural_basis
[0].matrix_space():
742 # Thanks for nothing! Matrix spaces aren't vector
743 # spaces in Sage, so we have to figure out its
744 # natural-basis coordinates ourselves.
745 V
= VectorSpace(elt
.base_ring(), elt
.nrows()**2)
746 W
= V
.span( _mat2vec(s
) for s
in natural_basis
)
747 coords
= W
.coordinates(_mat2vec(elt
))
748 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, coords
)
750 def __pow__(self
, n
):
752 Return ``self`` raised to the power ``n``.
754 Jordan algebras are always power-associative; see for
755 example Faraut and Koranyi, Proposition II.1.2 (ii).
759 We have to override this because our superclass uses row vectors
760 instead of column vectors! We, on the other hand, assume column
765 sage: set_random_seed()
766 sage: x = random_eja().random_element()
767 sage: x.operator()(x) == (x^2)
770 A few examples of power-associativity::
772 sage: set_random_seed()
773 sage: x = random_eja().random_element()
774 sage: x*(x*x)*(x*x) == x^5
776 sage: (x*x)*(x*x*x) == x^5
779 We also know that powers operator-commute (Koecher, Chapter
782 sage: set_random_seed()
783 sage: x = random_eja().random_element()
784 sage: m = ZZ.random_element(0,10)
785 sage: n = ZZ.random_element(0,10)
786 sage: Lxm = (x^m).operator()
787 sage: Lxn = (x^n).operator()
788 sage: Lxm*Lxn == Lxn*Lxm
793 return self
.parent().one()
797 return (self
.operator()**(n
-1))(self
)
800 def apply_univariate_polynomial(self
, p
):
802 Apply the univariate polynomial ``p`` to this element.
804 A priori, SageMath won't allow us to apply a univariate
805 polynomial to an element of an EJA, because we don't know
806 that EJAs are rings (they are usually not associative). Of
807 course, we know that EJAs are power-associative, so the
808 operation is ultimately kosher. This function sidesteps
809 the CAS to get the answer we want and expect.
813 sage: R = PolynomialRing(QQ, 't')
815 sage: p = t^4 - t^3 + 5*t - 2
816 sage: J = RealCartesianProductEJA(5)
817 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
822 We should always get back an element of the algebra::
824 sage: set_random_seed()
825 sage: p = PolynomialRing(QQ, 't').random_element()
826 sage: J = random_eja()
827 sage: x = J.random_element()
828 sage: x.apply_univariate_polynomial(p) in J
832 if len(p
.variables()) > 1:
833 raise ValueError("not a univariate polynomial")
836 # Convert the coeficcients to the parent's base ring,
837 # because a priori they might live in an (unnecessarily)
838 # larger ring for which P.sum() would fail below.
839 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
840 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
843 def characteristic_polynomial(self
):
845 Return the characteristic polynomial of this element.
849 The rank of `R^3` is three, and the minimal polynomial of
850 the identity element is `(t-1)` from which it follows that
851 the characteristic polynomial should be `(t-1)^3`::
853 sage: J = RealCartesianProductEJA(3)
854 sage: J.one().characteristic_polynomial()
855 t^3 - 3*t^2 + 3*t - 1
857 Likewise, the characteristic of the zero element in the
858 rank-three algebra `R^{n}` should be `t^{3}`::
860 sage: J = RealCartesianProductEJA(3)
861 sage: J.zero().characteristic_polynomial()
864 The characteristic polynomial of an element should evaluate
865 to zero on that element::
867 sage: set_random_seed()
868 sage: x = RealCartesianProductEJA(3).random_element()
869 sage: p = x.characteristic_polynomial()
870 sage: x.apply_univariate_polynomial(p)
874 p
= self
.parent().characteristic_polynomial()
875 return p(*self
.vector())
878 def inner_product(self
, other
):
880 Return the parent algebra's inner product of myself and ``other``.
884 The inner product in the Jordan spin algebra is the usual
885 inner product on `R^n` (this example only works because the
886 basis for the Jordan algebra is the standard basis in `R^n`)::
888 sage: J = JordanSpinEJA(3)
889 sage: x = vector(QQ,[1,2,3])
890 sage: y = vector(QQ,[4,5,6])
891 sage: x.inner_product(y)
893 sage: J(x).inner_product(J(y))
896 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
897 multiplication is the usual matrix multiplication in `S^n`,
898 so the inner product of the identity matrix with itself
901 sage: J = RealSymmetricEJA(3)
902 sage: J.one().inner_product(J.one())
905 Likewise, the inner product on `C^n` is `<X,Y> =
906 Re(trace(X*Y))`, where we must necessarily take the real
907 part because the product of Hermitian matrices may not be
910 sage: J = ComplexHermitianEJA(3)
911 sage: J.one().inner_product(J.one())
914 Ditto for the quaternions::
916 sage: J = QuaternionHermitianEJA(3)
917 sage: J.one().inner_product(J.one())
922 Ensure that we can always compute an inner product, and that
923 it gives us back a real number::
925 sage: set_random_seed()
926 sage: J = random_eja()
927 sage: x = J.random_element()
928 sage: y = J.random_element()
929 sage: x.inner_product(y) in RR
935 raise TypeError("'other' must live in the same algebra")
937 return P
.inner_product(self
, other
)
940 def operator_commutes_with(self
, other
):
942 Return whether or not this element operator-commutes
947 The definition of a Jordan algebra says that any element
948 operator-commutes with its square::
950 sage: set_random_seed()
951 sage: x = random_eja().random_element()
952 sage: x.operator_commutes_with(x^2)
957 Test Lemma 1 from Chapter III of Koecher::
959 sage: set_random_seed()
960 sage: J = random_eja()
961 sage: u = J.random_element()
962 sage: v = J.random_element()
963 sage: lhs = u.operator_commutes_with(u*v)
964 sage: rhs = v.operator_commutes_with(u^2)
968 Test the first polarization identity from my notes, Koecher Chapter
969 III, or from Baes (2.3)::
971 sage: set_random_seed()
972 sage: J = random_eja()
973 sage: x = J.random_element()
974 sage: y = J.random_element()
975 sage: Lx = x.operator()
976 sage: Ly = y.operator()
977 sage: Lxx = (x*x).operator()
978 sage: Lxy = (x*y).operator()
979 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
982 Test the second polarization identity from my notes or from
985 sage: set_random_seed()
986 sage: J = random_eja()
987 sage: x = J.random_element()
988 sage: y = J.random_element()
989 sage: z = J.random_element()
990 sage: Lx = x.operator()
991 sage: Ly = y.operator()
992 sage: Lz = z.operator()
993 sage: Lzy = (z*y).operator()
994 sage: Lxy = (x*y).operator()
995 sage: Lxz = (x*z).operator()
996 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
999 Test the third polarization identity from my notes or from
1002 sage: set_random_seed()
1003 sage: J = random_eja()
1004 sage: u = J.random_element()
1005 sage: y = J.random_element()
1006 sage: z = J.random_element()
1007 sage: Lu = u.operator()
1008 sage: Ly = y.operator()
1009 sage: Lz = z.operator()
1010 sage: Lzy = (z*y).operator()
1011 sage: Luy = (u*y).operator()
1012 sage: Luz = (u*z).operator()
1013 sage: Luyz = (u*(y*z)).operator()
1014 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
1015 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
1016 sage: bool(lhs == rhs)
1020 if not other
in self
.parent():
1021 raise TypeError("'other' must live in the same algebra")
1024 B
= other
.operator()
1030 Return my determinant, the product of my eigenvalues.
1034 sage: J = JordanSpinEJA(2)
1035 sage: e0,e1 = J.gens()
1036 sage: x = sum( J.gens() )
1042 sage: J = JordanSpinEJA(3)
1043 sage: e0,e1,e2 = J.gens()
1044 sage: x = sum( J.gens() )
1050 An element is invertible if and only if its determinant is
1053 sage: set_random_seed()
1054 sage: x = random_eja().random_element()
1055 sage: x.is_invertible() == (x.det() != 0)
1061 p
= P
._charpoly
_coeff
(0)
1062 # The _charpoly_coeff function already adds the factor of
1063 # -1 to ensure that _charpoly_coeff(0) is really what
1064 # appears in front of t^{0} in the charpoly. However,
1065 # we want (-1)^r times THAT for the determinant.
1066 return ((-1)**r
)*p(*self
.vector())
1071 Return the Jordan-multiplicative inverse of this element.
1075 We appeal to the quadratic representation as in Koecher's
1076 Theorem 12 in Chapter III, Section 5.
1080 The inverse in the spin factor algebra is given in Alizadeh's
1083 sage: set_random_seed()
1084 sage: n = ZZ.random_element(1,10)
1085 sage: J = JordanSpinEJA(n)
1086 sage: x = J.random_element()
1087 sage: while not x.is_invertible():
1088 ....: x = J.random_element()
1089 sage: x_vec = x.vector()
1091 sage: x_bar = x_vec[1:]
1092 sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
1093 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
1094 sage: x_inverse = coeff*inv_vec
1095 sage: x.inverse() == J(x_inverse)
1100 The identity element is its own inverse::
1102 sage: set_random_seed()
1103 sage: J = random_eja()
1104 sage: J.one().inverse() == J.one()
1107 If an element has an inverse, it acts like one::
1109 sage: set_random_seed()
1110 sage: J = random_eja()
1111 sage: x = J.random_element()
1112 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
1115 The inverse of the inverse is what we started with::
1117 sage: set_random_seed()
1118 sage: J = random_eja()
1119 sage: x = J.random_element()
1120 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
1123 The zero element is never invertible::
1125 sage: set_random_seed()
1126 sage: J = random_eja().zero().inverse()
1127 Traceback (most recent call last):
1129 ValueError: element is not invertible
1132 if not self
.is_invertible():
1133 raise ValueError("element is not invertible")
1135 return (~self
.quadratic_representation())(self
)
1138 def is_invertible(self
):
1140 Return whether or not this element is invertible.
1142 We can't use the superclass method because it relies on
1143 the algebra being associative.
1147 The usual way to do this is to check if the determinant is
1148 zero, but we need the characteristic polynomial for the
1149 determinant. The minimal polynomial is a lot easier to get,
1150 so we use Corollary 2 in Chapter V of Koecher to check
1151 whether or not the paren't algebra's zero element is a root
1152 of this element's minimal polynomial.
1156 The identity element is always invertible::
1158 sage: set_random_seed()
1159 sage: J = random_eja()
1160 sage: J.one().is_invertible()
1163 The zero element is never invertible::
1165 sage: set_random_seed()
1166 sage: J = random_eja()
1167 sage: J.zero().is_invertible()
1171 zero
= self
.parent().zero()
1172 p
= self
.minimal_polynomial()
1173 return not (p(zero
) == zero
)
1176 def is_nilpotent(self
):
1178 Return whether or not some power of this element is zero.
1180 The superclass method won't work unless we're in an
1181 associative algebra, and we aren't. However, we generate
1182 an assocoative subalgebra and we're nilpotent there if and
1183 only if we're nilpotent here (probably).
1187 The identity element is never nilpotent::
1189 sage: set_random_seed()
1190 sage: random_eja().one().is_nilpotent()
1193 The additive identity is always nilpotent::
1195 sage: set_random_seed()
1196 sage: random_eja().zero().is_nilpotent()
1200 # The element we're going to call "is_nilpotent()" on.
1201 # Either myself, interpreted as an element of a finite-
1202 # dimensional algebra, or an element of an associative
1206 if self
.parent().is_associative():
1207 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1209 V
= self
.span_of_powers()
1210 assoc_subalg
= self
.subalgebra_generated_by()
1211 # Mis-design warning: the basis used for span_of_powers()
1212 # and subalgebra_generated_by() must be the same, and in
1214 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1216 # Recursive call, but should work since elt lives in an
1217 # associative algebra.
1218 return elt
.is_nilpotent()
1221 def is_regular(self
):
1223 Return whether or not this is a regular element.
1227 The identity element always has degree one, but any element
1228 linearly-independent from it is regular::
1230 sage: J = JordanSpinEJA(5)
1231 sage: J.one().is_regular()
1233 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
1234 sage: for x in J.gens():
1235 ....: (J.one() + x).is_regular()
1243 return self
.degree() == self
.parent().rank()
1248 Compute the degree of this element the straightforward way
1249 according to the definition; by appending powers to a list
1250 and figuring out its dimension (that is, whether or not
1251 they're linearly dependent).
1255 sage: J = JordanSpinEJA(4)
1256 sage: J.one().degree()
1258 sage: e0,e1,e2,e3 = J.gens()
1259 sage: (e0 - e1).degree()
1262 In the spin factor algebra (of rank two), all elements that
1263 aren't multiples of the identity are regular::
1265 sage: set_random_seed()
1266 sage: n = ZZ.random_element(1,10)
1267 sage: J = JordanSpinEJA(n)
1268 sage: x = J.random_element()
1269 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
1273 return self
.span_of_powers().dimension()
1276 def left_matrix(self
):
1278 Our parent class defines ``left_matrix`` and ``matrix``
1279 methods whose names are misleading. We don't want them.
1281 raise NotImplementedError("use operator().matrix() instead")
1283 matrix
= left_matrix
1286 def minimal_polynomial(self
):
1288 Return the minimal polynomial of this element,
1289 as a function of the variable `t`.
1293 We restrict ourselves to the associative subalgebra
1294 generated by this element, and then return the minimal
1295 polynomial of this element's operator matrix (in that
1296 subalgebra). This works by Baes Proposition 2.3.16.
1300 The minimal polynomial of the identity and zero elements are
1303 sage: set_random_seed()
1304 sage: J = random_eja()
1305 sage: J.one().minimal_polynomial()
1307 sage: J.zero().minimal_polynomial()
1310 The degree of an element is (by one definition) the degree
1311 of its minimal polynomial::
1313 sage: set_random_seed()
1314 sage: x = random_eja().random_element()
1315 sage: x.degree() == x.minimal_polynomial().degree()
1318 The minimal polynomial and the characteristic polynomial coincide
1319 and are known (see Alizadeh, Example 11.11) for all elements of
1320 the spin factor algebra that aren't scalar multiples of the
1323 sage: set_random_seed()
1324 sage: n = ZZ.random_element(2,10)
1325 sage: J = JordanSpinEJA(n)
1326 sage: y = J.random_element()
1327 sage: while y == y.coefficient(0)*J.one():
1328 ....: y = J.random_element()
1329 sage: y0 = y.vector()[0]
1330 sage: y_bar = y.vector()[1:]
1331 sage: actual = y.minimal_polynomial()
1332 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1333 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1334 sage: bool(actual == expected)
1337 The minimal polynomial should always kill its element::
1339 sage: set_random_seed()
1340 sage: x = random_eja().random_element()
1341 sage: p = x.minimal_polynomial()
1342 sage: x.apply_univariate_polynomial(p)
1346 V
= self
.span_of_powers()
1347 assoc_subalg
= self
.subalgebra_generated_by()
1348 # Mis-design warning: the basis used for span_of_powers()
1349 # and subalgebra_generated_by() must be the same, and in
1351 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1353 # We get back a symbolic polynomial in 'x' but want a real
1354 # polynomial in 't'.
1355 p_of_x
= elt
.operator().matrix().minimal_polynomial()
1356 return p_of_x
.change_variable_name('t')
1359 def natural_representation(self
):
1361 Return a more-natural representation of this element.
1363 Every finite-dimensional Euclidean Jordan Algebra is a
1364 direct sum of five simple algebras, four of which comprise
1365 Hermitian matrices. This method returns the original
1366 "natural" representation of this element as a Hermitian
1367 matrix, if it has one. If not, you get the usual representation.
1371 sage: J = ComplexHermitianEJA(3)
1374 sage: J.one().natural_representation()
1384 sage: J = QuaternionHermitianEJA(3)
1387 sage: J.one().natural_representation()
1388 [1 0 0 0 0 0 0 0 0 0 0 0]
1389 [0 1 0 0 0 0 0 0 0 0 0 0]
1390 [0 0 1 0 0 0 0 0 0 0 0 0]
1391 [0 0 0 1 0 0 0 0 0 0 0 0]
1392 [0 0 0 0 1 0 0 0 0 0 0 0]
1393 [0 0 0 0 0 1 0 0 0 0 0 0]
1394 [0 0 0 0 0 0 1 0 0 0 0 0]
1395 [0 0 0 0 0 0 0 1 0 0 0 0]
1396 [0 0 0 0 0 0 0 0 1 0 0 0]
1397 [0 0 0 0 0 0 0 0 0 1 0 0]
1398 [0 0 0 0 0 0 0 0 0 0 1 0]
1399 [0 0 0 0 0 0 0 0 0 0 0 1]
1402 B
= self
.parent().natural_basis()
1403 W
= B
[0].matrix_space()
1404 return W
.linear_combination(zip(self
.vector(), B
))
1409 Return the left-multiplication-by-this-element
1410 operator on the ambient algebra.
1414 sage: set_random_seed()
1415 sage: J = random_eja()
1416 sage: x = J.random_element()
1417 sage: y = J.random_element()
1418 sage: x.operator()(y) == x*y
1420 sage: y.operator()(x) == x*y
1425 fda_elt
= FiniteDimensionalAlgebraElement(P
, self
)
1426 return FiniteDimensionalEuclideanJordanAlgebraOperator(
1429 fda_elt
.matrix().transpose() )
1432 def quadratic_representation(self
, other
=None):
1434 Return the quadratic representation of this element.
1438 The explicit form in the spin factor algebra is given by
1439 Alizadeh's Example 11.12::
1441 sage: set_random_seed()
1442 sage: n = ZZ.random_element(1,10)
1443 sage: J = JordanSpinEJA(n)
1444 sage: x = J.random_element()
1445 sage: x_vec = x.vector()
1447 sage: x_bar = x_vec[1:]
1448 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1449 sage: B = 2*x0*x_bar.row()
1450 sage: C = 2*x0*x_bar.column()
1451 sage: D = identity_matrix(QQ, n-1)
1452 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1453 sage: D = D + 2*x_bar.tensor_product(x_bar)
1454 sage: Q = block_matrix(2,2,[A,B,C,D])
1455 sage: Q == x.quadratic_representation().matrix()
1458 Test all of the properties from Theorem 11.2 in Alizadeh::
1460 sage: set_random_seed()
1461 sage: J = random_eja()
1462 sage: x = J.random_element()
1463 sage: y = J.random_element()
1464 sage: Lx = x.operator()
1465 sage: Lxx = (x*x).operator()
1466 sage: Qx = x.quadratic_representation()
1467 sage: Qy = y.quadratic_representation()
1468 sage: Qxy = x.quadratic_representation(y)
1469 sage: Qex = J.one().quadratic_representation(x)
1470 sage: n = ZZ.random_element(10)
1471 sage: Qxn = (x^n).quadratic_representation()
1475 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1478 Property 2 (multiply on the right for :trac:`28272`):
1480 sage: alpha = QQ.random_element()
1481 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1486 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1489 sage: not x.is_invertible() or (
1492 ....: x.inverse().quadratic_representation() )
1495 sage: Qxy(J.one()) == x*y
1500 sage: not x.is_invertible() or (
1501 ....: x.quadratic_representation(x.inverse())*Qx
1502 ....: == Qx*x.quadratic_representation(x.inverse()) )
1505 sage: not x.is_invertible() or (
1506 ....: x.quadratic_representation(x.inverse())*Qx
1508 ....: 2*x.operator()*Qex - Qx )
1511 sage: 2*x.operator()*Qex - Qx == Lxx
1516 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1526 sage: not x.is_invertible() or (
1527 ....: Qx*x.inverse().operator() == Lx )
1532 sage: not x.operator_commutes_with(y) or (
1533 ....: Qx(y)^n == Qxn(y^n) )
1539 elif not other
in self
.parent():
1540 raise TypeError("'other' must live in the same algebra")
1543 M
= other
.operator()
1544 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1547 def span_of_powers(self
):
1549 Return the vector space spanned by successive powers of
1552 # The dimension of the subalgebra can't be greater than
1553 # the big algebra, so just put everything into a list
1554 # and let span() get rid of the excess.
1556 # We do the extra ambient_vector_space() in case we're messing
1557 # with polynomials and the direct parent is a module.
1558 V
= self
.parent().vector_space()
1559 return V
.span( (self
**d
).vector() for d
in xrange(V
.dimension()) )
1562 def subalgebra_generated_by(self
):
1564 Return the associative subalgebra of the parent EJA generated
1569 sage: set_random_seed()
1570 sage: x = random_eja().random_element()
1571 sage: x.subalgebra_generated_by().is_associative()
1574 Squaring in the subalgebra should work the same as in
1577 sage: set_random_seed()
1578 sage: x = random_eja().random_element()
1579 sage: u = x.subalgebra_generated_by().random_element()
1580 sage: u.operator()(u) == u^2
1584 # First get the subspace spanned by the powers of myself...
1585 V
= self
.span_of_powers()
1586 F
= self
.base_ring()
1588 # Now figure out the entries of the right-multiplication
1589 # matrix for the successive basis elements b0, b1,... of
1592 for b_right
in V
.basis():
1593 eja_b_right
= self
.parent()(b_right
)
1595 # The first row of the right-multiplication matrix by
1596 # b1 is what we get if we apply that matrix to b1. The
1597 # second row of the right multiplication matrix by b1
1598 # is what we get when we apply that matrix to b2...
1600 # IMPORTANT: this assumes that all vectors are COLUMN
1601 # vectors, unlike our superclass (which uses row vectors).
1602 for b_left
in V
.basis():
1603 eja_b_left
= self
.parent()(b_left
)
1604 # Multiply in the original EJA, but then get the
1605 # coordinates from the subalgebra in terms of its
1607 this_row
= V
.coordinates((eja_b_left
*eja_b_right
).vector())
1608 b_right_rows
.append(this_row
)
1609 b_right_matrix
= matrix(F
, b_right_rows
)
1610 mats
.append(b_right_matrix
)
1612 # It's an algebra of polynomials in one element, and EJAs
1613 # are power-associative.
1615 # TODO: choose generator names intelligently.
1616 return FiniteDimensionalEuclideanJordanAlgebra(F
, mats
, assume_associative
=True, names
='f')
1619 def subalgebra_idempotent(self
):
1621 Find an idempotent in the associative subalgebra I generate
1622 using Proposition 2.3.5 in Baes.
1626 sage: set_random_seed()
1627 sage: J = random_eja()
1628 sage: x = J.random_element()
1629 sage: while x.is_nilpotent():
1630 ....: x = J.random_element()
1631 sage: c = x.subalgebra_idempotent()
1636 if self
.is_nilpotent():
1637 raise ValueError("this only works with non-nilpotent elements!")
1639 V
= self
.span_of_powers()
1640 J
= self
.subalgebra_generated_by()
1641 # Mis-design warning: the basis used for span_of_powers()
1642 # and subalgebra_generated_by() must be the same, and in
1644 u
= J(V
.coordinates(self
.vector()))
1646 # The image of the matrix of left-u^m-multiplication
1647 # will be minimal for some natural number s...
1649 minimal_dim
= V
.dimension()
1650 for i
in xrange(1, V
.dimension()):
1651 this_dim
= (u
**i
).operator().matrix().image().dimension()
1652 if this_dim
< minimal_dim
:
1653 minimal_dim
= this_dim
1656 # Now minimal_matrix should correspond to the smallest
1657 # non-zero subspace in Baes's (or really, Koecher's)
1660 # However, we need to restrict the matrix to work on the
1661 # subspace... or do we? Can't we just solve, knowing that
1662 # A(c) = u^(s+1) should have a solution in the big space,
1665 # Beware, solve_right() means that we're using COLUMN vectors.
1666 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1668 A
= u_next
.operator().matrix()
1669 c_coordinates
= A
.solve_right(u_next
.vector())
1671 # Now c_coordinates is the idempotent we want, but it's in
1672 # the coordinate system of the subalgebra.
1674 # We need the basis for J, but as elements of the parent algebra.
1676 basis
= [self
.parent(v
) for v
in V
.basis()]
1677 return self
.parent().linear_combination(zip(c_coordinates
, basis
))
1682 Return my trace, the sum of my eigenvalues.
1686 sage: J = JordanSpinEJA(3)
1687 sage: x = sum(J.gens())
1693 sage: J = RealCartesianProductEJA(5)
1694 sage: J.one().trace()
1699 The trace of an element is a real number::
1701 sage: set_random_seed()
1702 sage: J = random_eja()
1703 sage: J.random_element().trace() in J.base_ring()
1709 p
= P
._charpoly
_coeff
(r
-1)
1710 # The _charpoly_coeff function already adds the factor of
1711 # -1 to ensure that _charpoly_coeff(r-1) is really what
1712 # appears in front of t^{r-1} in the charpoly. However,
1713 # we want the negative of THAT for the trace.
1714 return -p(*self
.vector())
1717 def trace_inner_product(self
, other
):
1719 Return the trace inner product of myself and ``other``.
1723 The trace inner product is commutative::
1725 sage: set_random_seed()
1726 sage: J = random_eja()
1727 sage: x = J.random_element(); y = J.random_element()
1728 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1731 The trace inner product is bilinear::
1733 sage: set_random_seed()
1734 sage: J = random_eja()
1735 sage: x = J.random_element()
1736 sage: y = J.random_element()
1737 sage: z = J.random_element()
1738 sage: a = QQ.random_element();
1739 sage: actual = (a*(x+z)).trace_inner_product(y)
1740 sage: expected = ( a*x.trace_inner_product(y) +
1741 ....: a*z.trace_inner_product(y) )
1742 sage: actual == expected
1744 sage: actual = x.trace_inner_product(a*(y+z))
1745 sage: expected = ( a*x.trace_inner_product(y) +
1746 ....: a*x.trace_inner_product(z) )
1747 sage: actual == expected
1750 The trace inner product satisfies the compatibility
1751 condition in the definition of a Euclidean Jordan algebra::
1753 sage: set_random_seed()
1754 sage: J = random_eja()
1755 sage: x = J.random_element()
1756 sage: y = J.random_element()
1757 sage: z = J.random_element()
1758 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1762 if not other
in self
.parent():
1763 raise TypeError("'other' must live in the same algebra")
1765 return (self
*other
).trace()
1768 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1770 Return the Euclidean Jordan Algebra corresponding to the set
1771 `R^n` under the Hadamard product.
1773 Note: this is nothing more than the Cartesian product of ``n``
1774 copies of the spin algebra. Once Cartesian product algebras
1775 are implemented, this can go.
1779 This multiplication table can be verified by hand::
1781 sage: J = RealCartesianProductEJA(3)
1782 sage: e0,e1,e2 = J.gens()
1798 def __classcall_private__(cls
, n
, field
=QQ
):
1799 # The FiniteDimensionalAlgebra constructor takes a list of
1800 # matrices, the ith representing right multiplication by the ith
1801 # basis element in the vector space. So if e_1 = (1,0,0), then
1802 # right (Hadamard) multiplication of x by e_1 picks out the first
1803 # component of x; and likewise for the ith basis element e_i.
1804 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
1805 for i
in xrange(n
) ]
1807 fdeja
= super(RealCartesianProductEJA
, cls
)
1808 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=n
)
1810 def inner_product(self
, x
, y
):
1811 return _usual_ip(x
,y
)
1816 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1820 For now, we choose a random natural number ``n`` (greater than zero)
1821 and then give you back one of the following:
1823 * The cartesian product of the rational numbers ``n`` times; this is
1824 ``QQ^n`` with the Hadamard product.
1826 * The Jordan spin algebra on ``QQ^n``.
1828 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1831 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1832 in the space of ``2n``-by-``2n`` real symmetric matrices.
1834 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1835 in the space of ``4n``-by-``4n`` real symmetric matrices.
1837 Later this might be extended to return Cartesian products of the
1843 Euclidean Jordan algebra of degree...
1847 # The max_n component lets us choose different upper bounds on the
1848 # value "n" that gets passed to the constructor. This is needed
1849 # because e.g. R^{10} is reasonable to test, while the Hermitian
1850 # 10-by-10 quaternion matrices are not.
1851 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
1853 (RealSymmetricEJA
, 5),
1854 (ComplexHermitianEJA
, 4),
1855 (QuaternionHermitianEJA
, 3)])
1856 n
= ZZ
.random_element(1, max_n
)
1857 return constructor(n
, field
=QQ
)
1861 def _real_symmetric_basis(n
, field
=QQ
):
1863 Return a basis for the space of real symmetric n-by-n matrices.
1865 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1869 for j
in xrange(i
+1):
1870 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1874 # Beware, orthogonal but not normalized!
1875 Sij
= Eij
+ Eij
.transpose()
1880 def _complex_hermitian_basis(n
, field
=QQ
):
1882 Returns a basis for the space of complex Hermitian n-by-n matrices.
1886 sage: set_random_seed()
1887 sage: n = ZZ.random_element(1,5)
1888 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1892 F
= QuadraticField(-1, 'I')
1895 # This is like the symmetric case, but we need to be careful:
1897 # * We want conjugate-symmetry, not just symmetry.
1898 # * The diagonal will (as a result) be real.
1902 for j
in xrange(i
+1):
1903 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1905 Sij
= _embed_complex_matrix(Eij
)
1908 # Beware, orthogonal but not normalized! The second one
1909 # has a minus because it's conjugated.
1910 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
1912 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
1917 def _quaternion_hermitian_basis(n
, field
=QQ
):
1919 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1923 sage: set_random_seed()
1924 sage: n = ZZ.random_element(1,5)
1925 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1929 Q
= QuaternionAlgebra(QQ
,-1,-1)
1932 # This is like the symmetric case, but we need to be careful:
1934 # * We want conjugate-symmetry, not just symmetry.
1935 # * The diagonal will (as a result) be real.
1939 for j
in xrange(i
+1):
1940 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1942 Sij
= _embed_quaternion_matrix(Eij
)
1945 # Beware, orthogonal but not normalized! The second,
1946 # third, and fourth ones have a minus because they're
1948 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
1950 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
1952 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
1954 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
1960 return vector(m
.base_ring(), m
.list())
1963 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
1965 def _multiplication_table_from_matrix_basis(basis
):
1967 At least three of the five simple Euclidean Jordan algebras have the
1968 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1969 multiplication on the right is matrix multiplication. Given a basis
1970 for the underlying matrix space, this function returns a
1971 multiplication table (obtained by looping through the basis
1972 elements) for an algebra of those matrices. A reordered copy
1973 of the basis is also returned to work around the fact that
1974 the ``span()`` in this function will change the order of the basis
1975 from what we think it is, to... something else.
1977 # In S^2, for example, we nominally have four coordinates even
1978 # though the space is of dimension three only. The vector space V
1979 # is supposed to hold the entire long vector, and the subspace W
1980 # of V will be spanned by the vectors that arise from symmetric
1981 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1982 field
= basis
[0].base_ring()
1983 dimension
= basis
[0].nrows()
1985 V
= VectorSpace(field
, dimension
**2)
1986 W
= V
.span( _mat2vec(s
) for s
in basis
)
1988 # Taking the span above reorders our basis (thanks, jerk!) so we
1989 # need to put our "matrix basis" in the same order as the
1990 # (reordered) vector basis.
1991 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
1995 # Brute force the multiplication-by-s matrix by looping
1996 # through all elements of the basis and doing the computation
1997 # to find out what the corresponding row should be. BEWARE:
1998 # these multiplication tables won't be symmetric! It therefore
1999 # becomes REALLY IMPORTANT that the underlying algebra
2000 # constructor uses ROW vectors and not COLUMN vectors. That's
2001 # why we're computing rows here and not columns.
2004 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
2005 Q_rows
.append(W
.coordinates(this_row
))
2006 Q
= matrix(field
, W
.dimension(), Q_rows
)
2012 def _embed_complex_matrix(M
):
2014 Embed the n-by-n complex matrix ``M`` into the space of real
2015 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2016 bi` to the block matrix ``[[a,b],[-b,a]]``.
2020 sage: F = QuadraticField(-1,'i')
2021 sage: x1 = F(4 - 2*i)
2022 sage: x2 = F(1 + 2*i)
2025 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2026 sage: _embed_complex_matrix(M)
2035 Embedding is a homomorphism (isomorphism, in fact)::
2037 sage: set_random_seed()
2038 sage: n = ZZ.random_element(5)
2039 sage: F = QuadraticField(-1, 'i')
2040 sage: X = random_matrix(F, n)
2041 sage: Y = random_matrix(F, n)
2042 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
2043 sage: expected = _embed_complex_matrix(X*Y)
2044 sage: actual == expected
2050 raise ValueError("the matrix 'M' must be square")
2051 field
= M
.base_ring()
2056 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
2058 # We can drop the imaginaries here.
2059 return block_matrix(field
.base_ring(), n
, blocks
)
2062 def _unembed_complex_matrix(M
):
2064 The inverse of _embed_complex_matrix().
2068 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2069 ....: [-2, 1, -4, 3],
2070 ....: [ 9, 10, 11, 12],
2071 ....: [-10, 9, -12, 11] ])
2072 sage: _unembed_complex_matrix(A)
2074 [ 10*i + 9 12*i + 11]
2078 Unembedding is the inverse of embedding::
2080 sage: set_random_seed()
2081 sage: F = QuadraticField(-1, 'i')
2082 sage: M = random_matrix(F, 3)
2083 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
2089 raise ValueError("the matrix 'M' must be square")
2090 if not n
.mod(2).is_zero():
2091 raise ValueError("the matrix 'M' must be a complex embedding")
2093 F
= QuadraticField(-1, 'i')
2096 # Go top-left to bottom-right (reading order), converting every
2097 # 2-by-2 block we see to a single complex element.
2099 for k
in xrange(n
/2):
2100 for j
in xrange(n
/2):
2101 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
2102 if submat
[0,0] != submat
[1,1]:
2103 raise ValueError('bad on-diagonal submatrix')
2104 if submat
[0,1] != -submat
[1,0]:
2105 raise ValueError('bad off-diagonal submatrix')
2106 z
= submat
[0,0] + submat
[0,1]*i
2109 return matrix(F
, n
/2, elements
)
2112 def _embed_quaternion_matrix(M
):
2114 Embed the n-by-n quaternion matrix ``M`` into the space of real
2115 matrices of size 4n-by-4n by first sending each quaternion entry
2116 `z = a + bi + cj + dk` to the block-complex matrix
2117 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
2122 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2123 sage: i,j,k = Q.gens()
2124 sage: x = 1 + 2*i + 3*j + 4*k
2125 sage: M = matrix(Q, 1, [[x]])
2126 sage: _embed_quaternion_matrix(M)
2132 Embedding is a homomorphism (isomorphism, in fact)::
2134 sage: set_random_seed()
2135 sage: n = ZZ.random_element(5)
2136 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2137 sage: X = random_matrix(Q, n)
2138 sage: Y = random_matrix(Q, n)
2139 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
2140 sage: expected = _embed_quaternion_matrix(X*Y)
2141 sage: actual == expected
2145 quaternions
= M
.base_ring()
2148 raise ValueError("the matrix 'M' must be square")
2150 F
= QuadraticField(-1, 'i')
2155 t
= z
.coefficient_tuple()
2160 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2161 [-c
+ d
*i
, a
- b
*i
]])
2162 blocks
.append(_embed_complex_matrix(cplx_matrix
))
2164 # We should have real entries by now, so use the realest field
2165 # we've got for the return value.
2166 return block_matrix(quaternions
.base_ring(), n
, blocks
)
2169 def _unembed_quaternion_matrix(M
):
2171 The inverse of _embed_quaternion_matrix().
2175 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2176 ....: [-2, 1, -4, 3],
2177 ....: [-3, 4, 1, -2],
2178 ....: [-4, -3, 2, 1]])
2179 sage: _unembed_quaternion_matrix(M)
2180 [1 + 2*i + 3*j + 4*k]
2184 Unembedding is the inverse of embedding::
2186 sage: set_random_seed()
2187 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2188 sage: M = random_matrix(Q, 3)
2189 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
2195 raise ValueError("the matrix 'M' must be square")
2196 if not n
.mod(4).is_zero():
2197 raise ValueError("the matrix 'M' must be a complex embedding")
2199 Q
= QuaternionAlgebra(QQ
,-1,-1)
2202 # Go top-left to bottom-right (reading order), converting every
2203 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2206 for l
in xrange(n
/4):
2207 for m
in xrange(n
/4):
2208 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
2209 if submat
[0,0] != submat
[1,1].conjugate():
2210 raise ValueError('bad on-diagonal submatrix')
2211 if submat
[0,1] != -submat
[1,0].conjugate():
2212 raise ValueError('bad off-diagonal submatrix')
2213 z
= submat
[0,0].real() + submat
[0,0].imag()*i
2214 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
2217 return matrix(Q
, n
/4, elements
)
2220 # The usual inner product on R^n.
2222 return x
.vector().inner_product(y
.vector())
2224 # The inner product used for the real symmetric simple EJA.
2225 # We keep it as a separate function because e.g. the complex
2226 # algebra uses the same inner product, except divided by 2.
2227 def _matrix_ip(X
,Y
):
2228 X_mat
= X
.natural_representation()
2229 Y_mat
= Y
.natural_representation()
2230 return (X_mat
*Y_mat
).trace()
2233 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2235 The rank-n simple EJA consisting of real symmetric n-by-n
2236 matrices, the usual symmetric Jordan product, and the trace inner
2237 product. It has dimension `(n^2 + n)/2` over the reals.
2241 sage: J = RealSymmetricEJA(2)
2242 sage: e0, e1, e2 = J.gens()
2252 The degree of this algebra is `(n^2 + n) / 2`::
2254 sage: set_random_seed()
2255 sage: n = ZZ.random_element(1,5)
2256 sage: J = RealSymmetricEJA(n)
2257 sage: J.degree() == (n^2 + n)/2
2260 The Jordan multiplication is what we think it is::
2262 sage: set_random_seed()
2263 sage: n = ZZ.random_element(1,5)
2264 sage: J = RealSymmetricEJA(n)
2265 sage: x = J.random_element()
2266 sage: y = J.random_element()
2267 sage: actual = (x*y).natural_representation()
2268 sage: X = x.natural_representation()
2269 sage: Y = y.natural_representation()
2270 sage: expected = (X*Y + Y*X)/2
2271 sage: actual == expected
2273 sage: J(expected) == x*y
2278 def __classcall_private__(cls
, n
, field
=QQ
):
2279 S
= _real_symmetric_basis(n
, field
=field
)
2280 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2282 fdeja
= super(RealSymmetricEJA
, cls
)
2283 return fdeja
.__classcall
_private
__(cls
,
2289 def inner_product(self
, x
, y
):
2290 return _matrix_ip(x
,y
)
2293 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2295 The rank-n simple EJA consisting of complex Hermitian n-by-n
2296 matrices over the real numbers, the usual symmetric Jordan product,
2297 and the real-part-of-trace inner product. It has dimension `n^2` over
2302 The degree of this algebra is `n^2`::
2304 sage: set_random_seed()
2305 sage: n = ZZ.random_element(1,5)
2306 sage: J = ComplexHermitianEJA(n)
2307 sage: J.degree() == n^2
2310 The Jordan multiplication is what we think it is::
2312 sage: set_random_seed()
2313 sage: n = ZZ.random_element(1,5)
2314 sage: J = ComplexHermitianEJA(n)
2315 sage: x = J.random_element()
2316 sage: y = J.random_element()
2317 sage: actual = (x*y).natural_representation()
2318 sage: X = x.natural_representation()
2319 sage: Y = y.natural_representation()
2320 sage: expected = (X*Y + Y*X)/2
2321 sage: actual == expected
2323 sage: J(expected) == x*y
2328 def __classcall_private__(cls
, n
, field
=QQ
):
2329 S
= _complex_hermitian_basis(n
)
2330 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2332 fdeja
= super(ComplexHermitianEJA
, cls
)
2333 return fdeja
.__classcall
_private
__(cls
,
2339 def inner_product(self
, x
, y
):
2340 # Since a+bi on the diagonal is represented as
2345 # we'll double-count the "a" entries if we take the trace of
2347 return _matrix_ip(x
,y
)/2
2350 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2352 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2353 matrices, the usual symmetric Jordan product, and the
2354 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2359 The degree of this algebra is `n^2`::
2361 sage: set_random_seed()
2362 sage: n = ZZ.random_element(1,5)
2363 sage: J = QuaternionHermitianEJA(n)
2364 sage: J.degree() == 2*(n^2) - n
2367 The Jordan multiplication is what we think it is::
2369 sage: set_random_seed()
2370 sage: n = ZZ.random_element(1,5)
2371 sage: J = QuaternionHermitianEJA(n)
2372 sage: x = J.random_element()
2373 sage: y = J.random_element()
2374 sage: actual = (x*y).natural_representation()
2375 sage: X = x.natural_representation()
2376 sage: Y = y.natural_representation()
2377 sage: expected = (X*Y + Y*X)/2
2378 sage: actual == expected
2380 sage: J(expected) == x*y
2385 def __classcall_private__(cls
, n
, field
=QQ
):
2386 S
= _quaternion_hermitian_basis(n
)
2387 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2389 fdeja
= super(QuaternionHermitianEJA
, cls
)
2390 return fdeja
.__classcall
_private
__(cls
,
2396 def inner_product(self
, x
, y
):
2397 # Since a+bi+cj+dk on the diagonal is represented as
2399 # a + bi +cj + dk = [ a b c d]
2404 # we'll quadruple-count the "a" entries if we take the trace of
2406 return _matrix_ip(x
,y
)/4
2409 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2411 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2412 with the usual inner product and jordan product ``x*y =
2413 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2418 This multiplication table can be verified by hand::
2420 sage: J = JordanSpinEJA(4)
2421 sage: e0,e1,e2,e3 = J.gens()
2439 def __classcall_private__(cls
, n
, field
=QQ
):
2441 id_matrix
= identity_matrix(field
, n
)
2443 ei
= id_matrix
.column(i
)
2444 Qi
= zero_matrix(field
, n
)
2446 Qi
.set_column(0, ei
)
2447 Qi
+= diagonal_matrix(n
, [ei
[0]]*n
)
2448 # The addition of the diagonal matrix adds an extra ei[0] in the
2449 # upper-left corner of the matrix.
2450 Qi
[0,0] = Qi
[0,0] * ~
field(2)
2453 # The rank of the spin algebra is two, unless we're in a
2454 # one-dimensional ambient space (because the rank is bounded by
2455 # the ambient dimension).
2456 fdeja
= super(JordanSpinEJA
, cls
)
2457 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=min(n
,2))
2459 def inner_product(self
, x
, y
):
2460 return _usual_ip(x
,y
)