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.
32 parent
= Hom(domain_eja
, codomain_eja
, VectorSpaces(F
))
34 # The Map initializer will set our parent to a homset, which
35 # is explicitly NOT what we want, because these ain't algebra
37 super(FiniteDimensionalEuclideanJordanAlgebraOperator
,self
).__init
__(parent
)
39 # Keep a matrix around to do all of the real work. It would
40 # be nice if we could use a VectorSpaceMorphism instead, but
41 # those use row vectors that we don't want to accidentally
42 # expose to our users.
48 Allow this operator to be called only on elements of an EJA.
52 sage: J = JordanSpinEJA(3)
53 sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
54 sage: id = identity_matrix(J.base_ring(), J.dimension())
55 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
60 return self
.codomain()(self
.matrix()*x
.vector())
63 def _add_(self
, other
):
65 Add the ``other`` EJA operator to this one.
69 When we add two EJA operators, we get another one back::
71 sage: J = RealSymmetricEJA(2)
72 sage: id = identity_matrix(J.base_ring(), J.dimension())
73 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
74 sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
76 Linear operator between finite-dimensional Euclidean Jordan
77 algebras represented by the matrix:
81 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
82 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
84 If you try to add two identical vector space operators but on
85 different EJAs, that should blow up::
87 sage: J1 = RealSymmetricEJA(2)
88 sage: J2 = JordanSpinEJA(3)
89 sage: id = identity_matrix(QQ, 3)
90 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,J1,id)
91 sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,J2,id)
93 Traceback (most recent call last):
95 TypeError: unsupported operand parent(s) for +: ...
98 return FiniteDimensionalEuclideanJordanAlgebraOperator(
101 self
.matrix() + other
.matrix())
104 def _composition_(self
, other
, homset
):
106 Compose two EJA operators to get another one (and NOT a formal
107 composite object) back.
111 sage: J1 = JordanSpinEJA(3)
112 sage: J2 = RealCartesianProductEJA(2)
113 sage: J3 = RealSymmetricEJA(1)
114 sage: mat1 = matrix(QQ, [[1,2,3],
116 sage: mat2 = matrix(QQ, [[7,8]])
117 sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,
120 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,
124 Linear operator between finite-dimensional Euclidean Jordan
125 algebras represented by the matrix:
127 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
128 Codomain: Euclidean Jordan algebra of degree 1 over Rational Field
131 return FiniteDimensionalEuclideanJordanAlgebraOperator(
134 self
.matrix()*other
.matrix())
137 def __eq__(self
, other
):
138 if self
.domain() != other
.domain():
140 if self
.codomain() != other
.codomain():
142 if self
.matrix() != other
.matrix():
146 def __invert__(self
):
148 Invert this EJA operator.
152 sage: J = RealSymmetricEJA(2)
153 sage: id = identity_matrix(J.base_ring(), J.dimension())
154 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
156 Linear operator between finite-dimensional Euclidean Jordan
157 algebras represented by the matrix:
161 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
162 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
165 return FiniteDimensionalEuclideanJordanAlgebraOperator(
171 def __mul__(self
, other
):
173 Compose two EJA operators, or scale myself by an element of the
174 ambient vector space.
176 We need to override the real ``__mul__`` function to prevent the
177 coercion framework from throwing an error when it fails to convert
178 a base ring element into a morphism.
182 We can scale an operator on a rational algebra by a rational number::
184 sage: J = RealSymmetricEJA(2)
185 sage: e0,e1,e2 = J.gens()
186 sage: x = 2*e0 + 4*e1 + 16*e2
188 Linear operator between finite-dimensional Euclidean Jordan algebras
189 represented by the matrix:
193 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
194 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
195 sage: x.operator()*(1/2)
196 Linear operator between finite-dimensional Euclidean Jordan algebras
197 represented by the matrix:
201 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
202 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
205 if other
in self
.codomain().base_ring():
206 return FiniteDimensionalEuclideanJordanAlgebraOperator(
211 # This should eventually delegate to _composition_ after performing
212 # some sanity checks for us.
213 mor
= super(FiniteDimensionalEuclideanJordanAlgebraOperator
,self
)
214 return mor
.__mul
__(other
)
219 Negate this EJA operator.
223 sage: J = RealSymmetricEJA(2)
224 sage: id = identity_matrix(J.base_ring(), J.dimension())
225 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
227 Linear operator between finite-dimensional Euclidean Jordan
228 algebras represented by the matrix:
232 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
233 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
236 return FiniteDimensionalEuclideanJordanAlgebraOperator(
242 def __pow__(self
, n
):
244 Raise this EJA operator to the power ``n``.
248 Ensure that we get back another EJA operator that can be added,
249 subtracted, et cetera::
251 sage: J = RealSymmetricEJA(2)
252 sage: id = identity_matrix(J.base_ring(), J.dimension())
253 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
254 sage: f^0 + f^1 + f^2
255 Linear operator between finite-dimensional Euclidean Jordan
256 algebras represented by the matrix:
260 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
261 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
267 # Raising a vector space morphism to the zero power gives
268 # you back a special IdentityMorphism that is useless to us.
269 rows
= self
.codomain().dimension()
270 cols
= self
.domain().dimension()
271 mat
= matrix
.identity(self
.base_ring(), rows
, cols
)
273 mat
= self
.matrix()**n
275 return FiniteDimensionalEuclideanJordanAlgebraOperator(
284 A text representation of this linear operator on a Euclidean
289 sage: J = JordanSpinEJA(2)
290 sage: id = identity_matrix(J.base_ring(), J.dimension())
291 sage: FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
292 Linear operator between finite-dimensional Euclidean Jordan
293 algebras represented by the matrix:
296 Domain: Euclidean Jordan algebra of degree 2 over Rational Field
297 Codomain: Euclidean Jordan algebra of degree 2 over Rational Field
300 msg
= ("Linear operator between finite-dimensional Euclidean Jordan "
301 "algebras represented by the matrix:\n",
305 return ''.join(msg
).format(self
.matrix(),
310 def _sub_(self
, other
):
312 Subtract ``other`` from this EJA operator.
316 sage: J = RealSymmetricEJA(2)
317 sage: id = identity_matrix(J.base_ring(),J.dimension())
318 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
320 Linear operator between finite-dimensional Euclidean Jordan
321 algebras represented by the matrix:
325 Domain: Euclidean Jordan algebra of degree 3 over Rational Field
326 Codomain: Euclidean Jordan algebra of degree 3 over Rational Field
329 return (self
+ (-other
))
334 Return the matrix representation of this operator with respect
335 to the default bases of its (co)domain.
339 sage: J = RealSymmetricEJA(2)
340 sage: mat = matrix(J.base_ring(), J.dimension(), range(9))
341 sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,mat)
351 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra
):
353 def __classcall_private__(cls
,
357 assume_associative
=False,
362 mult_table
= [b
.base_extend(field
) for b
in mult_table
]
365 if not (is_Matrix(b
) and b
.dimensions() == (n
, n
)):
366 raise ValueError("input is not a multiplication table")
367 mult_table
= tuple(mult_table
)
369 cat
= FiniteDimensionalAlgebrasWithBasis(field
)
370 cat
.or_subcategory(category
)
371 if assume_associative
:
372 cat
= cat
.Associative()
374 names
= normalize_names(n
, names
)
376 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, cls
)
377 return fda
.__classcall
__(cls
,
380 assume_associative
=assume_associative
,
384 natural_basis
=natural_basis
)
391 assume_associative
=False,
398 By definition, Jordan multiplication commutes::
400 sage: set_random_seed()
401 sage: J = random_eja()
402 sage: x = J.random_element()
403 sage: y = J.random_element()
409 self
._natural
_basis
= natural_basis
410 self
._multiplication
_table
= mult_table
411 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
420 Return a string representation of ``self``.
422 fmt
= "Euclidean Jordan algebra of degree {} over {}"
423 return fmt
.format(self
.degree(), self
.base_ring())
426 def _a_regular_element(self
):
428 Guess a regular element. Needed to compute the basis for our
429 characteristic polynomial coefficients.
432 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
433 if not z
.is_regular():
434 raise ValueError("don't know a regular element")
439 def _charpoly_basis_space(self
):
441 Return the vector space spanned by the basis used in our
442 characteristic polynomial coefficients. This is used not only to
443 compute those coefficients, but also any time we need to
444 evaluate the coefficients (like when we compute the trace or
447 z
= self
._a
_regular
_element
()
448 V
= self
.vector_space()
449 V1
= V
.span_of_basis( (z
**k
).vector() for k
in range(self
.rank()) )
450 b
= (V1
.basis() + V1
.complement().basis())
451 return V
.span_of_basis(b
)
455 def _charpoly_coeff(self
, i
):
457 Return the coefficient polynomial "a_{i}" of this algebra's
458 general characteristic polynomial.
460 Having this be a separate cached method lets us compute and
461 store the trace/determinant (a_{r-1} and a_{0} respectively)
462 separate from the entire characteristic polynomial.
464 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
465 R
= A_of_x
.base_ring()
467 # Guaranteed by theory
470 # Danger: the in-place modification is done for performance
471 # reasons (reconstructing a matrix with huge polynomial
472 # entries is slow), but I don't know how cached_method works,
473 # so it's highly possible that we're modifying some global
474 # list variable by reference, here. In other words, you
475 # probably shouldn't call this method twice on the same
476 # algebra, at the same time, in two threads
477 Ai_orig
= A_of_x
.column(i
)
478 A_of_x
.set_column(i
,xr
)
479 numerator
= A_of_x
.det()
480 A_of_x
.set_column(i
,Ai_orig
)
482 # We're relying on the theory here to ensure that each a_i is
483 # indeed back in R, and the added negative signs are to make
484 # the whole charpoly expression sum to zero.
485 return R(-numerator
/detA
)
489 def _charpoly_matrix_system(self
):
491 Compute the matrix whose entries A_ij are polynomials in
492 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
493 corresponding to `x^r` and the determinent of the matrix A =
494 [A_ij]. In other words, all of the fixed (cachable) data needed
495 to compute the coefficients of the characteristic polynomial.
500 # Construct a new algebra over a multivariate polynomial ring...
501 names
= ['X' + str(i
) for i
in range(1,n
+1)]
502 R
= PolynomialRing(self
.base_ring(), names
)
503 J
= FiniteDimensionalEuclideanJordanAlgebra(R
,
504 self
._multiplication
_table
,
507 idmat
= identity_matrix(J
.base_ring(), n
)
509 W
= self
._charpoly
_basis
_space
()
510 W
= W
.change_ring(R
.fraction_field())
512 # Starting with the standard coordinates x = (X1,X2,...,Xn)
513 # and then converting the entries to W-coordinates allows us
514 # to pass in the standard coordinates to the charpoly and get
515 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
518 # W.coordinates(x^2) eval'd at (standard z-coords)
522 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
524 # We want the middle equivalent thing in our matrix, but use
525 # the first equivalent thing instead so that we can pass in
526 # standard coordinates.
527 x
= J(vector(R
, R
.gens()))
528 l1
= [column_matrix(W
.coordinates((x
**k
).vector())) for k
in range(r
)]
529 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
530 A_of_x
= block_matrix(R
, 1, n
, (l1
+ l2
))
531 xr
= W
.coordinates((x
**r
).vector())
532 return (A_of_x
, x
, xr
, A_of_x
.det())
536 def characteristic_polynomial(self
):
541 This implementation doesn't guarantee that the polynomial
542 denominator in the coefficients is not identically zero, so
543 theoretically it could crash. The way that this is handled
544 in e.g. Faraut and Koranyi is to use a basis that guarantees
545 the denominator is non-zero. But, doing so requires knowledge
546 of at least one regular element, and we don't even know how
547 to do that. The trade-off is that, if we use the standard basis,
548 the resulting polynomial will accept the "usual" coordinates. In
549 other words, we don't have to do a change of basis before e.g.
550 computing the trace or determinant.
554 The characteristic polynomial in the spin algebra is given in
555 Alizadeh, Example 11.11::
557 sage: J = JordanSpinEJA(3)
558 sage: p = J.characteristic_polynomial(); p
559 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
560 sage: xvec = J.one().vector()
568 # The list of coefficient polynomials a_1, a_2, ..., a_n.
569 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
571 # We go to a bit of trouble here to reorder the
572 # indeterminates, so that it's easier to evaluate the
573 # characteristic polynomial at x's coordinates and get back
574 # something in terms of t, which is what we want.
576 S
= PolynomialRing(self
.base_ring(),'t')
578 S
= PolynomialRing(S
, R
.variable_names())
581 # Note: all entries past the rth should be zero. The
582 # coefficient of the highest power (x^r) is 1, but it doesn't
583 # appear in the solution vector which contains coefficients
584 # for the other powers (to make them sum to x^r).
586 a
[r
] = 1 # corresponds to x^r
588 # When the rank is equal to the dimension, trying to
589 # assign a[r] goes out-of-bounds.
590 a
.append(1) # corresponds to x^r
592 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
595 def inner_product(self
, x
, y
):
597 The inner product associated with this Euclidean Jordan algebra.
599 Defaults to the trace inner product, but can be overridden by
600 subclasses if they are sure that the necessary properties are
605 The inner product must satisfy its axiom for this algebra to truly
606 be a Euclidean Jordan Algebra::
608 sage: set_random_seed()
609 sage: J = random_eja()
610 sage: x = J.random_element()
611 sage: y = J.random_element()
612 sage: z = J.random_element()
613 sage: (x*y).inner_product(z) == y.inner_product(x*z)
617 if (not x
in self
) or (not y
in self
):
618 raise TypeError("arguments must live in this algebra")
619 return x
.trace_inner_product(y
)
622 def natural_basis(self
):
624 Return a more-natural representation of this algebra's basis.
626 Every finite-dimensional Euclidean Jordan Algebra is a direct
627 sum of five simple algebras, four of which comprise Hermitian
628 matrices. This method returns the original "natural" basis
629 for our underlying vector space. (Typically, the natural basis
630 is used to construct the multiplication table in the first place.)
632 Note that this will always return a matrix. The standard basis
633 in `R^n` will be returned as `n`-by-`1` column matrices.
637 sage: J = RealSymmetricEJA(2)
640 sage: J.natural_basis()
648 sage: J = JordanSpinEJA(2)
651 sage: J.natural_basis()
658 if self
._natural
_basis
is None:
659 return tuple( b
.vector().column() for b
in self
.basis() )
661 return self
._natural
_basis
666 Return the rank of this EJA.
668 if self
._rank
is None:
669 raise ValueError("no rank specified at genesis")
673 def vector_space(self
):
675 Return the vector space that underlies this algebra.
679 sage: J = RealSymmetricEJA(2)
680 sage: J.vector_space()
681 Vector space of dimension 3 over Rational Field
684 return self
.zero().vector().parent().ambient_vector_space()
687 class Element(FiniteDimensionalAlgebraElement
):
689 An element of a Euclidean Jordan algebra.
694 Oh man, I should not be doing this. This hides the "disabled"
695 methods ``left_matrix`` and ``matrix`` from introspection;
696 in particular it removes them from tab-completion.
698 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
699 dir(self
.__class
__) )
702 def __init__(self
, A
, elt
=None):
706 The identity in `S^n` is converted to the identity in the EJA::
708 sage: J = RealSymmetricEJA(3)
709 sage: I = identity_matrix(QQ,3)
710 sage: J(I) == J.one()
713 This skew-symmetric matrix can't be represented in the EJA::
715 sage: J = RealSymmetricEJA(3)
716 sage: A = matrix(QQ,3, lambda i,j: i-j)
718 Traceback (most recent call last):
720 ArithmeticError: vector is not in free module
723 # Goal: if we're given a matrix, and if it lives in our
724 # parent algebra's "natural ambient space," convert it
725 # into an algebra element.
727 # The catch is, we make a recursive call after converting
728 # the given matrix into a vector that lives in the algebra.
729 # This we need to try the parent class initializer first,
730 # to avoid recursing forever if we're given something that
731 # already fits into the algebra, but also happens to live
732 # in the parent's "natural ambient space" (this happens with
735 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, elt
)
737 natural_basis
= A
.natural_basis()
738 if elt
in natural_basis
[0].matrix_space():
739 # Thanks for nothing! Matrix spaces aren't vector
740 # spaces in Sage, so we have to figure out its
741 # natural-basis coordinates ourselves.
742 V
= VectorSpace(elt
.base_ring(), elt
.nrows()**2)
743 W
= V
.span( _mat2vec(s
) for s
in natural_basis
)
744 coords
= W
.coordinates(_mat2vec(elt
))
745 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, coords
)
747 def __pow__(self
, n
):
749 Return ``self`` raised to the power ``n``.
751 Jordan algebras are always power-associative; see for
752 example Faraut and Koranyi, Proposition II.1.2 (ii).
756 We have to override this because our superclass uses row vectors
757 instead of column vectors! We, on the other hand, assume column
762 sage: set_random_seed()
763 sage: x = random_eja().random_element()
764 sage: x.operator_matrix()*x.vector() == (x^2).vector()
767 A few examples of power-associativity::
769 sage: set_random_seed()
770 sage: x = random_eja().random_element()
771 sage: x*(x*x)*(x*x) == x^5
773 sage: (x*x)*(x*x*x) == x^5
776 We also know that powers operator-commute (Koecher, Chapter
779 sage: set_random_seed()
780 sage: x = random_eja().random_element()
781 sage: m = ZZ.random_element(0,10)
782 sage: n = ZZ.random_element(0,10)
783 sage: Lxm = (x^m).operator_matrix()
784 sage: Lxn = (x^n).operator_matrix()
785 sage: Lxm*Lxn == Lxn*Lxm
795 return A( (self
.operator_matrix()**(n
-1))*self
.vector() )
798 def apply_univariate_polynomial(self
, p
):
800 Apply the univariate polynomial ``p`` to this element.
802 A priori, SageMath won't allow us to apply a univariate
803 polynomial to an element of an EJA, because we don't know
804 that EJAs are rings (they are usually not associative). Of
805 course, we know that EJAs are power-associative, so the
806 operation is ultimately kosher. This function sidesteps
807 the CAS to get the answer we want and expect.
811 sage: R = PolynomialRing(QQ, 't')
813 sage: p = t^4 - t^3 + 5*t - 2
814 sage: J = RealCartesianProductEJA(5)
815 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
820 We should always get back an element of the algebra::
822 sage: set_random_seed()
823 sage: p = PolynomialRing(QQ, 't').random_element()
824 sage: J = random_eja()
825 sage: x = J.random_element()
826 sage: x.apply_univariate_polynomial(p) in J
830 if len(p
.variables()) > 1:
831 raise ValueError("not a univariate polynomial")
834 # Convert the coeficcients to the parent's base ring,
835 # because a priori they might live in an (unnecessarily)
836 # larger ring for which P.sum() would fail below.
837 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
838 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
841 def characteristic_polynomial(self
):
843 Return the characteristic polynomial of this element.
847 The rank of `R^3` is three, and the minimal polynomial of
848 the identity element is `(t-1)` from which it follows that
849 the characteristic polynomial should be `(t-1)^3`::
851 sage: J = RealCartesianProductEJA(3)
852 sage: J.one().characteristic_polynomial()
853 t^3 - 3*t^2 + 3*t - 1
855 Likewise, the characteristic of the zero element in the
856 rank-three algebra `R^{n}` should be `t^{3}`::
858 sage: J = RealCartesianProductEJA(3)
859 sage: J.zero().characteristic_polynomial()
862 The characteristic polynomial of an element should evaluate
863 to zero on that element::
865 sage: set_random_seed()
866 sage: x = RealCartesianProductEJA(3).random_element()
867 sage: p = x.characteristic_polynomial()
868 sage: x.apply_univariate_polynomial(p)
872 p
= self
.parent().characteristic_polynomial()
873 return p(*self
.vector())
876 def inner_product(self
, other
):
878 Return the parent algebra's inner product of myself and ``other``.
882 The inner product in the Jordan spin algebra is the usual
883 inner product on `R^n` (this example only works because the
884 basis for the Jordan algebra is the standard basis in `R^n`)::
886 sage: J = JordanSpinEJA(3)
887 sage: x = vector(QQ,[1,2,3])
888 sage: y = vector(QQ,[4,5,6])
889 sage: x.inner_product(y)
891 sage: J(x).inner_product(J(y))
894 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
895 multiplication is the usual matrix multiplication in `S^n`,
896 so the inner product of the identity matrix with itself
899 sage: J = RealSymmetricEJA(3)
900 sage: J.one().inner_product(J.one())
903 Likewise, the inner product on `C^n` is `<X,Y> =
904 Re(trace(X*Y))`, where we must necessarily take the real
905 part because the product of Hermitian matrices may not be
908 sage: J = ComplexHermitianEJA(3)
909 sage: J.one().inner_product(J.one())
912 Ditto for the quaternions::
914 sage: J = QuaternionHermitianEJA(3)
915 sage: J.one().inner_product(J.one())
920 Ensure that we can always compute an inner product, and that
921 it gives us back a real number::
923 sage: set_random_seed()
924 sage: J = random_eja()
925 sage: x = J.random_element()
926 sage: y = J.random_element()
927 sage: x.inner_product(y) in RR
933 raise TypeError("'other' must live in the same algebra")
935 return P
.inner_product(self
, other
)
938 def operator_commutes_with(self
, other
):
940 Return whether or not this element operator-commutes
945 The definition of a Jordan algebra says that any element
946 operator-commutes with its square::
948 sage: set_random_seed()
949 sage: x = random_eja().random_element()
950 sage: x.operator_commutes_with(x^2)
955 Test Lemma 1 from Chapter III of Koecher::
957 sage: set_random_seed()
958 sage: J = random_eja()
959 sage: u = J.random_element()
960 sage: v = J.random_element()
961 sage: lhs = u.operator_commutes_with(u*v)
962 sage: rhs = v.operator_commutes_with(u^2)
967 if not other
in self
.parent():
968 raise TypeError("'other' must live in the same algebra")
970 A
= self
.operator_matrix()
971 B
= other
.operator_matrix()
977 Return my determinant, the product of my eigenvalues.
981 sage: J = JordanSpinEJA(2)
982 sage: e0,e1 = J.gens()
983 sage: x = sum( J.gens() )
989 sage: J = JordanSpinEJA(3)
990 sage: e0,e1,e2 = J.gens()
991 sage: x = sum( J.gens() )
997 An element is invertible if and only if its determinant is
1000 sage: set_random_seed()
1001 sage: x = random_eja().random_element()
1002 sage: x.is_invertible() == (x.det() != 0)
1008 p
= P
._charpoly
_coeff
(0)
1009 # The _charpoly_coeff function already adds the factor of
1010 # -1 to ensure that _charpoly_coeff(0) is really what
1011 # appears in front of t^{0} in the charpoly. However,
1012 # we want (-1)^r times THAT for the determinant.
1013 return ((-1)**r
)*p(*self
.vector())
1018 Return the Jordan-multiplicative inverse of this element.
1022 We appeal to the quadratic representation as in Koecher's
1023 Theorem 12 in Chapter III, Section 5.
1027 The inverse in the spin factor algebra is given in Alizadeh's
1030 sage: set_random_seed()
1031 sage: n = ZZ.random_element(1,10)
1032 sage: J = JordanSpinEJA(n)
1033 sage: x = J.random_element()
1034 sage: while not x.is_invertible():
1035 ....: x = J.random_element()
1036 sage: x_vec = x.vector()
1038 sage: x_bar = x_vec[1:]
1039 sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
1040 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
1041 sage: x_inverse = coeff*inv_vec
1042 sage: x.inverse() == J(x_inverse)
1047 The identity element is its own inverse::
1049 sage: set_random_seed()
1050 sage: J = random_eja()
1051 sage: J.one().inverse() == J.one()
1054 If an element has an inverse, it acts like one::
1056 sage: set_random_seed()
1057 sage: J = random_eja()
1058 sage: x = J.random_element()
1059 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
1062 The inverse of the inverse is what we started with::
1064 sage: set_random_seed()
1065 sage: J = random_eja()
1066 sage: x = J.random_element()
1067 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
1070 The zero element is never invertible::
1072 sage: set_random_seed()
1073 sage: J = random_eja().zero().inverse()
1074 Traceback (most recent call last):
1076 ValueError: element is not invertible
1079 if not self
.is_invertible():
1080 raise ValueError("element is not invertible")
1082 return (~self
.quadratic_representation())(self
)
1085 def is_invertible(self
):
1087 Return whether or not this element is invertible.
1089 We can't use the superclass method because it relies on
1090 the algebra being associative.
1094 The usual way to do this is to check if the determinant is
1095 zero, but we need the characteristic polynomial for the
1096 determinant. The minimal polynomial is a lot easier to get,
1097 so we use Corollary 2 in Chapter V of Koecher to check
1098 whether or not the paren't algebra's zero element is a root
1099 of this element's minimal polynomial.
1103 The identity element is always invertible::
1105 sage: set_random_seed()
1106 sage: J = random_eja()
1107 sage: J.one().is_invertible()
1110 The zero element is never invertible::
1112 sage: set_random_seed()
1113 sage: J = random_eja()
1114 sage: J.zero().is_invertible()
1118 zero
= self
.parent().zero()
1119 p
= self
.minimal_polynomial()
1120 return not (p(zero
) == zero
)
1123 def is_nilpotent(self
):
1125 Return whether or not some power of this element is zero.
1127 The superclass method won't work unless we're in an
1128 associative algebra, and we aren't. However, we generate
1129 an assocoative subalgebra and we're nilpotent there if and
1130 only if we're nilpotent here (probably).
1134 The identity element is never nilpotent::
1136 sage: set_random_seed()
1137 sage: random_eja().one().is_nilpotent()
1140 The additive identity is always nilpotent::
1142 sage: set_random_seed()
1143 sage: random_eja().zero().is_nilpotent()
1147 # The element we're going to call "is_nilpotent()" on.
1148 # Either myself, interpreted as an element of a finite-
1149 # dimensional algebra, or an element of an associative
1153 if self
.parent().is_associative():
1154 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1156 V
= self
.span_of_powers()
1157 assoc_subalg
= self
.subalgebra_generated_by()
1158 # Mis-design warning: the basis used for span_of_powers()
1159 # and subalgebra_generated_by() must be the same, and in
1161 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1163 # Recursive call, but should work since elt lives in an
1164 # associative algebra.
1165 return elt
.is_nilpotent()
1168 def is_regular(self
):
1170 Return whether or not this is a regular element.
1174 The identity element always has degree one, but any element
1175 linearly-independent from it is regular::
1177 sage: J = JordanSpinEJA(5)
1178 sage: J.one().is_regular()
1180 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
1181 sage: for x in J.gens():
1182 ....: (J.one() + x).is_regular()
1190 return self
.degree() == self
.parent().rank()
1195 Compute the degree of this element the straightforward way
1196 according to the definition; by appending powers to a list
1197 and figuring out its dimension (that is, whether or not
1198 they're linearly dependent).
1202 sage: J = JordanSpinEJA(4)
1203 sage: J.one().degree()
1205 sage: e0,e1,e2,e3 = J.gens()
1206 sage: (e0 - e1).degree()
1209 In the spin factor algebra (of rank two), all elements that
1210 aren't multiples of the identity are regular::
1212 sage: set_random_seed()
1213 sage: n = ZZ.random_element(1,10)
1214 sage: J = JordanSpinEJA(n)
1215 sage: x = J.random_element()
1216 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
1220 return self
.span_of_powers().dimension()
1223 def left_matrix(self
):
1225 Our parent class defines ``left_matrix`` and ``matrix``
1226 methods whose names are misleading. We don't want them.
1228 raise NotImplementedError("use operator_matrix() instead")
1230 matrix
= left_matrix
1233 def minimal_polynomial(self
):
1235 Return the minimal polynomial of this element,
1236 as a function of the variable `t`.
1240 We restrict ourselves to the associative subalgebra
1241 generated by this element, and then return the minimal
1242 polynomial of this element's operator matrix (in that
1243 subalgebra). This works by Baes Proposition 2.3.16.
1247 The minimal polynomial of the identity and zero elements are
1250 sage: set_random_seed()
1251 sage: J = random_eja()
1252 sage: J.one().minimal_polynomial()
1254 sage: J.zero().minimal_polynomial()
1257 The degree of an element is (by one definition) the degree
1258 of its minimal polynomial::
1260 sage: set_random_seed()
1261 sage: x = random_eja().random_element()
1262 sage: x.degree() == x.minimal_polynomial().degree()
1265 The minimal polynomial and the characteristic polynomial coincide
1266 and are known (see Alizadeh, Example 11.11) for all elements of
1267 the spin factor algebra that aren't scalar multiples of the
1270 sage: set_random_seed()
1271 sage: n = ZZ.random_element(2,10)
1272 sage: J = JordanSpinEJA(n)
1273 sage: y = J.random_element()
1274 sage: while y == y.coefficient(0)*J.one():
1275 ....: y = J.random_element()
1276 sage: y0 = y.vector()[0]
1277 sage: y_bar = y.vector()[1:]
1278 sage: actual = y.minimal_polynomial()
1279 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1280 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1281 sage: bool(actual == expected)
1284 The minimal polynomial should always kill its element::
1286 sage: set_random_seed()
1287 sage: x = random_eja().random_element()
1288 sage: p = x.minimal_polynomial()
1289 sage: x.apply_univariate_polynomial(p)
1293 V
= self
.span_of_powers()
1294 assoc_subalg
= self
.subalgebra_generated_by()
1295 # Mis-design warning: the basis used for span_of_powers()
1296 # and subalgebra_generated_by() must be the same, and in
1298 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1300 # We get back a symbolic polynomial in 'x' but want a real
1301 # polynomial in 't'.
1302 p_of_x
= elt
.operator_matrix().minimal_polynomial()
1303 return p_of_x
.change_variable_name('t')
1306 def natural_representation(self
):
1308 Return a more-natural representation of this element.
1310 Every finite-dimensional Euclidean Jordan Algebra is a
1311 direct sum of five simple algebras, four of which comprise
1312 Hermitian matrices. This method returns the original
1313 "natural" representation of this element as a Hermitian
1314 matrix, if it has one. If not, you get the usual representation.
1318 sage: J = ComplexHermitianEJA(3)
1321 sage: J.one().natural_representation()
1331 sage: J = QuaternionHermitianEJA(3)
1334 sage: J.one().natural_representation()
1335 [1 0 0 0 0 0 0 0 0 0 0 0]
1336 [0 1 0 0 0 0 0 0 0 0 0 0]
1337 [0 0 1 0 0 0 0 0 0 0 0 0]
1338 [0 0 0 1 0 0 0 0 0 0 0 0]
1339 [0 0 0 0 1 0 0 0 0 0 0 0]
1340 [0 0 0 0 0 1 0 0 0 0 0 0]
1341 [0 0 0 0 0 0 1 0 0 0 0 0]
1342 [0 0 0 0 0 0 0 1 0 0 0 0]
1343 [0 0 0 0 0 0 0 0 1 0 0 0]
1344 [0 0 0 0 0 0 0 0 0 1 0 0]
1345 [0 0 0 0 0 0 0 0 0 0 1 0]
1346 [0 0 0 0 0 0 0 0 0 0 0 1]
1349 B
= self
.parent().natural_basis()
1350 W
= B
[0].matrix_space()
1351 return W
.linear_combination(zip(self
.vector(), B
))
1356 Return the left-multiplication-by-this-element
1357 operator on the ambient algebra.
1361 sage: set_random_seed()
1362 sage: J = random_eja()
1363 sage: x = J.random_element()
1364 sage: y = J.random_element()
1365 sage: x.operator()(y) == x*y
1367 sage: y.operator()(x) == x*y
1372 return FiniteDimensionalEuclideanJordanAlgebraOperator(
1374 self
.operator_matrix() )
1378 def operator_matrix(self
):
1380 Return the matrix that represents left- (or right-)
1381 multiplication by this element in the parent algebra.
1383 We implement this ourselves to work around the fact that
1384 our parent class represents everything with row vectors.
1388 Test the first polarization identity from my notes, Koecher Chapter
1389 III, or from Baes (2.3)::
1391 sage: set_random_seed()
1392 sage: J = random_eja()
1393 sage: x = J.random_element()
1394 sage: y = J.random_element()
1395 sage: Lx = x.operator_matrix()
1396 sage: Ly = y.operator_matrix()
1397 sage: Lxx = (x*x).operator_matrix()
1398 sage: Lxy = (x*y).operator_matrix()
1399 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
1402 Test the second polarization identity from my notes or from
1405 sage: set_random_seed()
1406 sage: J = random_eja()
1407 sage: x = J.random_element()
1408 sage: y = J.random_element()
1409 sage: z = J.random_element()
1410 sage: Lx = x.operator_matrix()
1411 sage: Ly = y.operator_matrix()
1412 sage: Lz = z.operator_matrix()
1413 sage: Lzy = (z*y).operator_matrix()
1414 sage: Lxy = (x*y).operator_matrix()
1415 sage: Lxz = (x*z).operator_matrix()
1416 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
1419 Test the third polarization identity from my notes or from
1422 sage: set_random_seed()
1423 sage: J = random_eja()
1424 sage: u = J.random_element()
1425 sage: y = J.random_element()
1426 sage: z = J.random_element()
1427 sage: Lu = u.operator_matrix()
1428 sage: Ly = y.operator_matrix()
1429 sage: Lz = z.operator_matrix()
1430 sage: Lzy = (z*y).operator_matrix()
1431 sage: Luy = (u*y).operator_matrix()
1432 sage: Luz = (u*z).operator_matrix()
1433 sage: Luyz = (u*(y*z)).operator_matrix()
1434 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
1435 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
1436 sage: bool(lhs == rhs)
1439 Ensure that our operator's ``matrix`` method agrees with
1440 this implementation::
1442 sage: set_random_seed()
1443 sage: J = random_eja()
1444 sage: x = J.random_element()
1445 sage: x.operator().matrix() == x.operator_matrix()
1449 fda_elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1450 return fda_elt
.matrix().transpose()
1453 def quadratic_representation(self
, other
=None):
1455 Return the quadratic representation of this element.
1459 The explicit form in the spin factor algebra is given by
1460 Alizadeh's Example 11.12::
1462 sage: set_random_seed()
1463 sage: n = ZZ.random_element(1,10)
1464 sage: J = JordanSpinEJA(n)
1465 sage: x = J.random_element()
1466 sage: x_vec = x.vector()
1468 sage: x_bar = x_vec[1:]
1469 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1470 sage: B = 2*x0*x_bar.row()
1471 sage: C = 2*x0*x_bar.column()
1472 sage: D = identity_matrix(QQ, n-1)
1473 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1474 sage: D = D + 2*x_bar.tensor_product(x_bar)
1475 sage: Q = block_matrix(2,2,[A,B,C,D])
1476 sage: Q == x.quadratic_representation().matrix()
1479 Test all of the properties from Theorem 11.2 in Alizadeh::
1481 sage: set_random_seed()
1482 sage: J = random_eja()
1483 sage: x = J.random_element()
1484 sage: y = J.random_element()
1485 sage: Lx = x.operator()
1486 sage: Lxx = (x*x).operator()
1487 sage: Qx = x.quadratic_representation()
1488 sage: Qy = y.quadratic_representation()
1489 sage: Qxy = x.quadratic_representation(y)
1490 sage: Qex = J.one().quadratic_representation(x)
1491 sage: n = ZZ.random_element(10)
1492 sage: Qxn = (x^n).quadratic_representation()
1496 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1499 Property 2 (multiply on the right for :trac:`28272`):
1501 sage: alpha = QQ.random_element()
1502 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1507 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1510 sage: not x.is_invertible() or (
1513 ....: x.inverse().quadratic_representation() )
1516 sage: Qxy(J.one()) == x*y
1521 sage: not x.is_invertible() or (
1522 ....: x.quadratic_representation(x.inverse())*Qx
1523 ....: == Qx*x.quadratic_representation(x.inverse()) )
1526 sage: not x.is_invertible() or (
1527 ....: x.quadratic_representation(x.inverse())*Qx
1529 ....: 2*x.operator()*Qex - Qx )
1532 sage: 2*x.operator()*Qex - Qx == Lxx
1537 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1547 sage: not x.is_invertible() or (
1548 ....: Qx*x.inverse().operator() == Lx )
1553 sage: not x.operator_commutes_with(y) or (
1554 ....: Qx(y)^n == Qxn(y^n) )
1560 elif not other
in self
.parent():
1561 raise TypeError("'other' must live in the same algebra")
1564 M
= other
.operator()
1565 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1568 def span_of_powers(self
):
1570 Return the vector space spanned by successive powers of
1573 # The dimension of the subalgebra can't be greater than
1574 # the big algebra, so just put everything into a list
1575 # and let span() get rid of the excess.
1577 # We do the extra ambient_vector_space() in case we're messing
1578 # with polynomials and the direct parent is a module.
1579 V
= self
.parent().vector_space()
1580 return V
.span( (self
**d
).vector() for d
in xrange(V
.dimension()) )
1583 def subalgebra_generated_by(self
):
1585 Return the associative subalgebra of the parent EJA generated
1590 sage: set_random_seed()
1591 sage: x = random_eja().random_element()
1592 sage: x.subalgebra_generated_by().is_associative()
1595 Squaring in the subalgebra should be the same thing as
1596 squaring in the superalgebra::
1598 sage: set_random_seed()
1599 sage: x = random_eja().random_element()
1600 sage: u = x.subalgebra_generated_by().random_element()
1601 sage: u.operator_matrix()*u.vector() == (u**2).vector()
1605 # First get the subspace spanned by the powers of myself...
1606 V
= self
.span_of_powers()
1607 F
= self
.base_ring()
1609 # Now figure out the entries of the right-multiplication
1610 # matrix for the successive basis elements b0, b1,... of
1613 for b_right
in V
.basis():
1614 eja_b_right
= self
.parent()(b_right
)
1616 # The first row of the right-multiplication matrix by
1617 # b1 is what we get if we apply that matrix to b1. The
1618 # second row of the right multiplication matrix by b1
1619 # is what we get when we apply that matrix to b2...
1621 # IMPORTANT: this assumes that all vectors are COLUMN
1622 # vectors, unlike our superclass (which uses row vectors).
1623 for b_left
in V
.basis():
1624 eja_b_left
= self
.parent()(b_left
)
1625 # Multiply in the original EJA, but then get the
1626 # coordinates from the subalgebra in terms of its
1628 this_row
= V
.coordinates((eja_b_left
*eja_b_right
).vector())
1629 b_right_rows
.append(this_row
)
1630 b_right_matrix
= matrix(F
, b_right_rows
)
1631 mats
.append(b_right_matrix
)
1633 # It's an algebra of polynomials in one element, and EJAs
1634 # are power-associative.
1636 # TODO: choose generator names intelligently.
1637 return FiniteDimensionalEuclideanJordanAlgebra(F
, mats
, assume_associative
=True, names
='f')
1640 def subalgebra_idempotent(self
):
1642 Find an idempotent in the associative subalgebra I generate
1643 using Proposition 2.3.5 in Baes.
1647 sage: set_random_seed()
1648 sage: J = random_eja()
1649 sage: x = J.random_element()
1650 sage: while x.is_nilpotent():
1651 ....: x = J.random_element()
1652 sage: c = x.subalgebra_idempotent()
1657 if self
.is_nilpotent():
1658 raise ValueError("this only works with non-nilpotent elements!")
1660 V
= self
.span_of_powers()
1661 J
= self
.subalgebra_generated_by()
1662 # Mis-design warning: the basis used for span_of_powers()
1663 # and subalgebra_generated_by() must be the same, and in
1665 u
= J(V
.coordinates(self
.vector()))
1667 # The image of the matrix of left-u^m-multiplication
1668 # will be minimal for some natural number s...
1670 minimal_dim
= V
.dimension()
1671 for i
in xrange(1, V
.dimension()):
1672 this_dim
= (u
**i
).operator_matrix().image().dimension()
1673 if this_dim
< minimal_dim
:
1674 minimal_dim
= this_dim
1677 # Now minimal_matrix should correspond to the smallest
1678 # non-zero subspace in Baes's (or really, Koecher's)
1681 # However, we need to restrict the matrix to work on the
1682 # subspace... or do we? Can't we just solve, knowing that
1683 # A(c) = u^(s+1) should have a solution in the big space,
1686 # Beware, solve_right() means that we're using COLUMN vectors.
1687 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1689 A
= u_next
.operator_matrix()
1690 c_coordinates
= A
.solve_right(u_next
.vector())
1692 # Now c_coordinates is the idempotent we want, but it's in
1693 # the coordinate system of the subalgebra.
1695 # We need the basis for J, but as elements of the parent algebra.
1697 basis
= [self
.parent(v
) for v
in V
.basis()]
1698 return self
.parent().linear_combination(zip(c_coordinates
, basis
))
1703 Return my trace, the sum of my eigenvalues.
1707 sage: J = JordanSpinEJA(3)
1708 sage: x = sum(J.gens())
1714 sage: J = RealCartesianProductEJA(5)
1715 sage: J.one().trace()
1720 The trace of an element is a real number::
1722 sage: set_random_seed()
1723 sage: J = random_eja()
1724 sage: J.random_element().trace() in J.base_ring()
1730 p
= P
._charpoly
_coeff
(r
-1)
1731 # The _charpoly_coeff function already adds the factor of
1732 # -1 to ensure that _charpoly_coeff(r-1) is really what
1733 # appears in front of t^{r-1} in the charpoly. However,
1734 # we want the negative of THAT for the trace.
1735 return -p(*self
.vector())
1738 def trace_inner_product(self
, other
):
1740 Return the trace inner product of myself and ``other``.
1744 The trace inner product is commutative::
1746 sage: set_random_seed()
1747 sage: J = random_eja()
1748 sage: x = J.random_element(); y = J.random_element()
1749 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1752 The trace inner product is bilinear::
1754 sage: set_random_seed()
1755 sage: J = random_eja()
1756 sage: x = J.random_element()
1757 sage: y = J.random_element()
1758 sage: z = J.random_element()
1759 sage: a = QQ.random_element();
1760 sage: actual = (a*(x+z)).trace_inner_product(y)
1761 sage: expected = ( a*x.trace_inner_product(y) +
1762 ....: a*z.trace_inner_product(y) )
1763 sage: actual == expected
1765 sage: actual = x.trace_inner_product(a*(y+z))
1766 sage: expected = ( a*x.trace_inner_product(y) +
1767 ....: a*x.trace_inner_product(z) )
1768 sage: actual == expected
1771 The trace inner product satisfies the compatibility
1772 condition in the definition of a Euclidean Jordan algebra::
1774 sage: set_random_seed()
1775 sage: J = random_eja()
1776 sage: x = J.random_element()
1777 sage: y = J.random_element()
1778 sage: z = J.random_element()
1779 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1783 if not other
in self
.parent():
1784 raise TypeError("'other' must live in the same algebra")
1786 return (self
*other
).trace()
1789 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1791 Return the Euclidean Jordan Algebra corresponding to the set
1792 `R^n` under the Hadamard product.
1794 Note: this is nothing more than the Cartesian product of ``n``
1795 copies of the spin algebra. Once Cartesian product algebras
1796 are implemented, this can go.
1800 This multiplication table can be verified by hand::
1802 sage: J = RealCartesianProductEJA(3)
1803 sage: e0,e1,e2 = J.gens()
1819 def __classcall_private__(cls
, n
, field
=QQ
):
1820 # The FiniteDimensionalAlgebra constructor takes a list of
1821 # matrices, the ith representing right multiplication by the ith
1822 # basis element in the vector space. So if e_1 = (1,0,0), then
1823 # right (Hadamard) multiplication of x by e_1 picks out the first
1824 # component of x; and likewise for the ith basis element e_i.
1825 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
1826 for i
in xrange(n
) ]
1828 fdeja
= super(RealCartesianProductEJA
, cls
)
1829 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=n
)
1831 def inner_product(self
, x
, y
):
1832 return _usual_ip(x
,y
)
1837 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1841 For now, we choose a random natural number ``n`` (greater than zero)
1842 and then give you back one of the following:
1844 * The cartesian product of the rational numbers ``n`` times; this is
1845 ``QQ^n`` with the Hadamard product.
1847 * The Jordan spin algebra on ``QQ^n``.
1849 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1852 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1853 in the space of ``2n``-by-``2n`` real symmetric matrices.
1855 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1856 in the space of ``4n``-by-``4n`` real symmetric matrices.
1858 Later this might be extended to return Cartesian products of the
1864 Euclidean Jordan algebra of degree...
1868 # The max_n component lets us choose different upper bounds on the
1869 # value "n" that gets passed to the constructor. This is needed
1870 # because e.g. R^{10} is reasonable to test, while the Hermitian
1871 # 10-by-10 quaternion matrices are not.
1872 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
1874 (RealSymmetricEJA
, 5),
1875 (ComplexHermitianEJA
, 4),
1876 (QuaternionHermitianEJA
, 3)])
1877 n
= ZZ
.random_element(1, max_n
)
1878 return constructor(n
, field
=QQ
)
1882 def _real_symmetric_basis(n
, field
=QQ
):
1884 Return a basis for the space of real symmetric n-by-n matrices.
1886 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1890 for j
in xrange(i
+1):
1891 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1895 # Beware, orthogonal but not normalized!
1896 Sij
= Eij
+ Eij
.transpose()
1901 def _complex_hermitian_basis(n
, field
=QQ
):
1903 Returns a basis for the space of complex Hermitian n-by-n matrices.
1907 sage: set_random_seed()
1908 sage: n = ZZ.random_element(1,5)
1909 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1913 F
= QuadraticField(-1, 'I')
1916 # This is like the symmetric case, but we need to be careful:
1918 # * We want conjugate-symmetry, not just symmetry.
1919 # * The diagonal will (as a result) be real.
1923 for j
in xrange(i
+1):
1924 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1926 Sij
= _embed_complex_matrix(Eij
)
1929 # Beware, orthogonal but not normalized! The second one
1930 # has a minus because it's conjugated.
1931 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
1933 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
1938 def _quaternion_hermitian_basis(n
, field
=QQ
):
1940 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1944 sage: set_random_seed()
1945 sage: n = ZZ.random_element(1,5)
1946 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1950 Q
= QuaternionAlgebra(QQ
,-1,-1)
1953 # This is like the symmetric case, but we need to be careful:
1955 # * We want conjugate-symmetry, not just symmetry.
1956 # * The diagonal will (as a result) be real.
1960 for j
in xrange(i
+1):
1961 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1963 Sij
= _embed_quaternion_matrix(Eij
)
1966 # Beware, orthogonal but not normalized! The second,
1967 # third, and fourth ones have a minus because they're
1969 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
1971 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
1973 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
1975 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
1981 return vector(m
.base_ring(), m
.list())
1984 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
1986 def _multiplication_table_from_matrix_basis(basis
):
1988 At least three of the five simple Euclidean Jordan algebras have the
1989 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1990 multiplication on the right is matrix multiplication. Given a basis
1991 for the underlying matrix space, this function returns a
1992 multiplication table (obtained by looping through the basis
1993 elements) for an algebra of those matrices. A reordered copy
1994 of the basis is also returned to work around the fact that
1995 the ``span()`` in this function will change the order of the basis
1996 from what we think it is, to... something else.
1998 # In S^2, for example, we nominally have four coordinates even
1999 # though the space is of dimension three only. The vector space V
2000 # is supposed to hold the entire long vector, and the subspace W
2001 # of V will be spanned by the vectors that arise from symmetric
2002 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
2003 field
= basis
[0].base_ring()
2004 dimension
= basis
[0].nrows()
2006 V
= VectorSpace(field
, dimension
**2)
2007 W
= V
.span( _mat2vec(s
) for s
in basis
)
2009 # Taking the span above reorders our basis (thanks, jerk!) so we
2010 # need to put our "matrix basis" in the same order as the
2011 # (reordered) vector basis.
2012 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
2016 # Brute force the multiplication-by-s matrix by looping
2017 # through all elements of the basis and doing the computation
2018 # to find out what the corresponding row should be. BEWARE:
2019 # these multiplication tables won't be symmetric! It therefore
2020 # becomes REALLY IMPORTANT that the underlying algebra
2021 # constructor uses ROW vectors and not COLUMN vectors. That's
2022 # why we're computing rows here and not columns.
2025 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
2026 Q_rows
.append(W
.coordinates(this_row
))
2027 Q
= matrix(field
, W
.dimension(), Q_rows
)
2033 def _embed_complex_matrix(M
):
2035 Embed the n-by-n complex matrix ``M`` into the space of real
2036 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2037 bi` to the block matrix ``[[a,b],[-b,a]]``.
2041 sage: F = QuadraticField(-1,'i')
2042 sage: x1 = F(4 - 2*i)
2043 sage: x2 = F(1 + 2*i)
2046 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2047 sage: _embed_complex_matrix(M)
2056 Embedding is a homomorphism (isomorphism, in fact)::
2058 sage: set_random_seed()
2059 sage: n = ZZ.random_element(5)
2060 sage: F = QuadraticField(-1, 'i')
2061 sage: X = random_matrix(F, n)
2062 sage: Y = random_matrix(F, n)
2063 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
2064 sage: expected = _embed_complex_matrix(X*Y)
2065 sage: actual == expected
2071 raise ValueError("the matrix 'M' must be square")
2072 field
= M
.base_ring()
2077 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
2079 # We can drop the imaginaries here.
2080 return block_matrix(field
.base_ring(), n
, blocks
)
2083 def _unembed_complex_matrix(M
):
2085 The inverse of _embed_complex_matrix().
2089 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2090 ....: [-2, 1, -4, 3],
2091 ....: [ 9, 10, 11, 12],
2092 ....: [-10, 9, -12, 11] ])
2093 sage: _unembed_complex_matrix(A)
2095 [ 10*i + 9 12*i + 11]
2099 Unembedding is the inverse of embedding::
2101 sage: set_random_seed()
2102 sage: F = QuadraticField(-1, 'i')
2103 sage: M = random_matrix(F, 3)
2104 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
2110 raise ValueError("the matrix 'M' must be square")
2111 if not n
.mod(2).is_zero():
2112 raise ValueError("the matrix 'M' must be a complex embedding")
2114 F
= QuadraticField(-1, 'i')
2117 # Go top-left to bottom-right (reading order), converting every
2118 # 2-by-2 block we see to a single complex element.
2120 for k
in xrange(n
/2):
2121 for j
in xrange(n
/2):
2122 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
2123 if submat
[0,0] != submat
[1,1]:
2124 raise ValueError('bad on-diagonal submatrix')
2125 if submat
[0,1] != -submat
[1,0]:
2126 raise ValueError('bad off-diagonal submatrix')
2127 z
= submat
[0,0] + submat
[0,1]*i
2130 return matrix(F
, n
/2, elements
)
2133 def _embed_quaternion_matrix(M
):
2135 Embed the n-by-n quaternion matrix ``M`` into the space of real
2136 matrices of size 4n-by-4n by first sending each quaternion entry
2137 `z = a + bi + cj + dk` to the block-complex matrix
2138 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
2143 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2144 sage: i,j,k = Q.gens()
2145 sage: x = 1 + 2*i + 3*j + 4*k
2146 sage: M = matrix(Q, 1, [[x]])
2147 sage: _embed_quaternion_matrix(M)
2153 Embedding is a homomorphism (isomorphism, in fact)::
2155 sage: set_random_seed()
2156 sage: n = ZZ.random_element(5)
2157 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2158 sage: X = random_matrix(Q, n)
2159 sage: Y = random_matrix(Q, n)
2160 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
2161 sage: expected = _embed_quaternion_matrix(X*Y)
2162 sage: actual == expected
2166 quaternions
= M
.base_ring()
2169 raise ValueError("the matrix 'M' must be square")
2171 F
= QuadraticField(-1, 'i')
2176 t
= z
.coefficient_tuple()
2181 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2182 [-c
+ d
*i
, a
- b
*i
]])
2183 blocks
.append(_embed_complex_matrix(cplx_matrix
))
2185 # We should have real entries by now, so use the realest field
2186 # we've got for the return value.
2187 return block_matrix(quaternions
.base_ring(), n
, blocks
)
2190 def _unembed_quaternion_matrix(M
):
2192 The inverse of _embed_quaternion_matrix().
2196 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2197 ....: [-2, 1, -4, 3],
2198 ....: [-3, 4, 1, -2],
2199 ....: [-4, -3, 2, 1]])
2200 sage: _unembed_quaternion_matrix(M)
2201 [1 + 2*i + 3*j + 4*k]
2205 Unembedding is the inverse of embedding::
2207 sage: set_random_seed()
2208 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2209 sage: M = random_matrix(Q, 3)
2210 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
2216 raise ValueError("the matrix 'M' must be square")
2217 if not n
.mod(4).is_zero():
2218 raise ValueError("the matrix 'M' must be a complex embedding")
2220 Q
= QuaternionAlgebra(QQ
,-1,-1)
2223 # Go top-left to bottom-right (reading order), converting every
2224 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2227 for l
in xrange(n
/4):
2228 for m
in xrange(n
/4):
2229 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
2230 if submat
[0,0] != submat
[1,1].conjugate():
2231 raise ValueError('bad on-diagonal submatrix')
2232 if submat
[0,1] != -submat
[1,0].conjugate():
2233 raise ValueError('bad off-diagonal submatrix')
2234 z
= submat
[0,0].real() + submat
[0,0].imag()*i
2235 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
2238 return matrix(Q
, n
/4, elements
)
2241 # The usual inner product on R^n.
2243 return x
.vector().inner_product(y
.vector())
2245 # The inner product used for the real symmetric simple EJA.
2246 # We keep it as a separate function because e.g. the complex
2247 # algebra uses the same inner product, except divided by 2.
2248 def _matrix_ip(X
,Y
):
2249 X_mat
= X
.natural_representation()
2250 Y_mat
= Y
.natural_representation()
2251 return (X_mat
*Y_mat
).trace()
2254 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2256 The rank-n simple EJA consisting of real symmetric n-by-n
2257 matrices, the usual symmetric Jordan product, and the trace inner
2258 product. It has dimension `(n^2 + n)/2` over the reals.
2262 sage: J = RealSymmetricEJA(2)
2263 sage: e0, e1, e2 = J.gens()
2273 The degree of this algebra is `(n^2 + n) / 2`::
2275 sage: set_random_seed()
2276 sage: n = ZZ.random_element(1,5)
2277 sage: J = RealSymmetricEJA(n)
2278 sage: J.degree() == (n^2 + n)/2
2281 The Jordan multiplication is what we think it is::
2283 sage: set_random_seed()
2284 sage: n = ZZ.random_element(1,5)
2285 sage: J = RealSymmetricEJA(n)
2286 sage: x = J.random_element()
2287 sage: y = J.random_element()
2288 sage: actual = (x*y).natural_representation()
2289 sage: X = x.natural_representation()
2290 sage: Y = y.natural_representation()
2291 sage: expected = (X*Y + Y*X)/2
2292 sage: actual == expected
2294 sage: J(expected) == x*y
2299 def __classcall_private__(cls
, n
, field
=QQ
):
2300 S
= _real_symmetric_basis(n
, field
=field
)
2301 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2303 fdeja
= super(RealSymmetricEJA
, cls
)
2304 return fdeja
.__classcall
_private
__(cls
,
2310 def inner_product(self
, x
, y
):
2311 return _matrix_ip(x
,y
)
2314 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2316 The rank-n simple EJA consisting of complex Hermitian n-by-n
2317 matrices over the real numbers, the usual symmetric Jordan product,
2318 and the real-part-of-trace inner product. It has dimension `n^2` over
2323 The degree of this algebra is `n^2`::
2325 sage: set_random_seed()
2326 sage: n = ZZ.random_element(1,5)
2327 sage: J = ComplexHermitianEJA(n)
2328 sage: J.degree() == n^2
2331 The Jordan multiplication is what we think it is::
2333 sage: set_random_seed()
2334 sage: n = ZZ.random_element(1,5)
2335 sage: J = ComplexHermitianEJA(n)
2336 sage: x = J.random_element()
2337 sage: y = J.random_element()
2338 sage: actual = (x*y).natural_representation()
2339 sage: X = x.natural_representation()
2340 sage: Y = y.natural_representation()
2341 sage: expected = (X*Y + Y*X)/2
2342 sage: actual == expected
2344 sage: J(expected) == x*y
2349 def __classcall_private__(cls
, n
, field
=QQ
):
2350 S
= _complex_hermitian_basis(n
)
2351 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2353 fdeja
= super(ComplexHermitianEJA
, cls
)
2354 return fdeja
.__classcall
_private
__(cls
,
2360 def inner_product(self
, x
, y
):
2361 # Since a+bi on the diagonal is represented as
2366 # we'll double-count the "a" entries if we take the trace of
2368 return _matrix_ip(x
,y
)/2
2371 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2373 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2374 matrices, the usual symmetric Jordan product, and the
2375 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2380 The degree of this algebra is `n^2`::
2382 sage: set_random_seed()
2383 sage: n = ZZ.random_element(1,5)
2384 sage: J = QuaternionHermitianEJA(n)
2385 sage: J.degree() == 2*(n^2) - n
2388 The Jordan multiplication is what we think it is::
2390 sage: set_random_seed()
2391 sage: n = ZZ.random_element(1,5)
2392 sage: J = QuaternionHermitianEJA(n)
2393 sage: x = J.random_element()
2394 sage: y = J.random_element()
2395 sage: actual = (x*y).natural_representation()
2396 sage: X = x.natural_representation()
2397 sage: Y = y.natural_representation()
2398 sage: expected = (X*Y + Y*X)/2
2399 sage: actual == expected
2401 sage: J(expected) == x*y
2406 def __classcall_private__(cls
, n
, field
=QQ
):
2407 S
= _quaternion_hermitian_basis(n
)
2408 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2410 fdeja
= super(QuaternionHermitianEJA
, cls
)
2411 return fdeja
.__classcall
_private
__(cls
,
2417 def inner_product(self
, x
, y
):
2418 # Since a+bi+cj+dk on the diagonal is represented as
2420 # a + bi +cj + dk = [ a b c d]
2425 # we'll quadruple-count the "a" entries if we take the trace of
2427 return _matrix_ip(x
,y
)/4
2430 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2432 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2433 with the usual inner product and jordan product ``x*y =
2434 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2439 This multiplication table can be verified by hand::
2441 sage: J = JordanSpinEJA(4)
2442 sage: e0,e1,e2,e3 = J.gens()
2460 def __classcall_private__(cls
, n
, field
=QQ
):
2462 id_matrix
= identity_matrix(field
, n
)
2464 ei
= id_matrix
.column(i
)
2465 Qi
= zero_matrix(field
, n
)
2467 Qi
.set_column(0, ei
)
2468 Qi
+= diagonal_matrix(n
, [ei
[0]]*n
)
2469 # The addition of the diagonal matrix adds an extra ei[0] in the
2470 # upper-left corner of the matrix.
2471 Qi
[0,0] = Qi
[0,0] * ~
field(2)
2474 # The rank of the spin algebra is two, unless we're in a
2475 # one-dimensional ambient space (because the rank is bounded by
2476 # the ambient dimension).
2477 fdeja
= super(JordanSpinEJA
, cls
)
2478 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=min(n
,2))
2480 def inner_product(self
, x
, y
):
2481 return _usual_ip(x
,y
)