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
.structure
.element
import is_Matrix
10 from sage
.structure
.category_object
import normalize_names
12 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra
import FiniteDimensionalAlgebra
13 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra_element
import FiniteDimensionalAlgebraElement
14 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra_morphism
import FiniteDimensionalAlgebraMorphism
17 class FiniteDimensionalEuclideanJordanAlgebraMorphism(FiniteDimensionalAlgebraMorphism
):
19 A linear map between two finite-dimensional EJAs.
21 This is a very thin wrapper around FiniteDimensionalAlgebraMorphism
22 that does only a few things:
24 1. Avoids the ``unitary`` and ``check`` arguments to the constructor
25 that will always be ``False``. This is necessary because these
26 are homomorphisms with respect to ADDITION, but the SageMath
27 machinery wants to check that they're homomorphisms with respect
28 to (Jordan) MULTIPLICATION. That obviously doesn't work.
30 2. Inputs and outputs the underlying matrix with respect to COLUMN
31 vectors, unlike the parent class.
33 3. Allows us to add, subtract, negate, multiply (compose), and
34 invert morphisms in the obvious way.
36 If this seems a bit heavyweight, it is. I would have been happy to
37 use a the ring morphism that underlies the finite-dimensional
38 algebra morphism, but they don't seem to be callable on elements of
39 our EJA, and you can't add/multiply/etc. them.
42 def _add_(self
, other
):
44 Add two EJA morphisms in the obvious way.
48 sage: J = RealSymmetricEJA(3)
51 sage: x.operator() + y.operator()
52 Morphism from Euclidean Jordan algebra of degree 6 over Rational
53 Field to Euclidean Jordan algebra of degree 6 over Rational Field
64 sage: set_random_seed()
65 sage: J = random_eja()
66 sage: x = J.random_element()
67 sage: y = J.random_element()
68 sage: (x.operator() + y.operator()) in J.Hom(J)
74 raise ValueError("summands must live in the same space")
76 return FiniteDimensionalEuclideanJordanAlgebraMorphism(
78 self
.matrix() + other
.matrix() )
81 def __init__(self
, parent
, f
):
82 FiniteDimensionalAlgebraMorphism
.__init
__(self
,
93 sage: J = RealSymmetricEJA(2)
94 sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
95 sage: x.is_invertible()
98 Morphism from Euclidean Jordan algebra of degree 3 over Rational
99 Field to Euclidean Jordan algebra of degree 3 over Rational Field
104 sage: x.operator_matrix().inverse()
111 sage: set_random_seed()
112 sage: J = random_eja()
113 sage: x = J.random_element()
114 sage: not x.is_invertible() or (
115 ....: (~x.operator()).matrix() == x.operator_matrix().inverse() )
120 if not A
.is_invertible():
121 raise ValueError("morphism is not invertible")
124 return FiniteDimensionalEuclideanJordanAlgebraMorphism(self
.parent(),
127 def _lmul_(self
, other
):
129 Compose two EJA morphisms using multiplicative notation.
133 sage: J = RealSymmetricEJA(3)
136 sage: x.operator() * y.operator()
137 Morphism from Euclidean Jordan algebra of degree 6 over Rational
138 Field to Euclidean Jordan algebra of degree 6 over Rational Field
149 sage: set_random_seed()
150 sage: J = random_eja()
151 sage: x = J.random_element()
152 sage: y = J.random_element()
153 sage: (x.operator() * y.operator()) in J.Hom(J)
157 if not other
.codomain() is self
.domain():
158 raise ValueError("(co)domains must agree for composition")
160 return FiniteDimensionalEuclideanJordanAlgebraMorphism(
162 self
.matrix()*other
.matrix() )
167 Negate this morphism.
171 sage: J = RealSymmetricEJA(2)
174 Morphism from Euclidean Jordan algebra of degree 3 over Rational
175 Field to Euclidean Jordan algebra of degree 3 over Rational Field
183 sage: set_random_seed()
184 sage: J = random_eja()
185 sage: x = J.random_element()
186 sage: -x.operator() in J.Hom(J)
190 return FiniteDimensionalEuclideanJordanAlgebraMorphism(
197 We override only the representation that is shown to the user,
198 because we want the matrix to be with respect to COLUMN vectors.
202 Ensure that we see the transpose of the underlying matrix object:
204 sage: J = RealSymmetricEJA(3)
205 sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
206 sage: L = x.operator()
208 Morphism from Euclidean Jordan algebra of degree 6 over Rational
209 Field to Euclidean Jordan algebra of degree 6 over Rational Field
226 return "Morphism from {} to {} given by matrix\n{}".format(
227 self
.domain(), self
.codomain(), self
.matrix())
230 def __sub__(self
, other
):
232 Subtract one morphism from another using addition and negation.
236 sage: J = RealSymmetricEJA(2)
237 sage: L1 = J.one().operator()
239 Morphism from Euclidean Jordan algebra of degree 3 over Rational
240 Field to Euclidean Jordan algebra of degree 3 over Rational
241 Field given by matrix
248 sage: set_random_seed()
249 sage: J = random_eja()
250 sage: x = J.random_element()
251 sage: y = J.random_element()
252 sage: x.operator() - y.operator() in J.Hom(J)
256 return self
+ (-other
)
261 Return the matrix of this morphism with respect to a left-action
264 return FiniteDimensionalAlgebraMorphism
.matrix(self
).transpose()
267 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra
):
269 def __classcall_private__(cls
,
273 assume_associative
=False,
278 mult_table
= [b
.base_extend(field
) for b
in mult_table
]
281 if not (is_Matrix(b
) and b
.dimensions() == (n
, n
)):
282 raise ValueError("input is not a multiplication table")
283 mult_table
= tuple(mult_table
)
285 cat
= FiniteDimensionalAlgebrasWithBasis(field
)
286 cat
.or_subcategory(category
)
287 if assume_associative
:
288 cat
= cat
.Associative()
290 names
= normalize_names(n
, names
)
292 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, cls
)
293 return fda
.__classcall
__(cls
,
296 assume_associative
=assume_associative
,
300 natural_basis
=natural_basis
)
307 assume_associative
=False,
314 By definition, Jordan multiplication commutes::
316 sage: set_random_seed()
317 sage: J = random_eja()
318 sage: x = J.random_element()
319 sage: y = J.random_element()
325 self
._natural
_basis
= natural_basis
326 self
._multiplication
_table
= mult_table
327 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
336 Return a string representation of ``self``.
338 fmt
= "Euclidean Jordan algebra of degree {} over {}"
339 return fmt
.format(self
.degree(), self
.base_ring())
342 def _a_regular_element(self
):
344 Guess a regular element. Needed to compute the basis for our
345 characteristic polynomial coefficients.
348 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
349 if not z
.is_regular():
350 raise ValueError("don't know a regular element")
355 def _charpoly_basis_space(self
):
357 Return the vector space spanned by the basis used in our
358 characteristic polynomial coefficients. This is used not only to
359 compute those coefficients, but also any time we need to
360 evaluate the coefficients (like when we compute the trace or
363 z
= self
._a
_regular
_element
()
364 V
= z
.vector().parent().ambient_vector_space()
365 V1
= V
.span_of_basis( (z
**k
).vector() for k
in range(self
.rank()) )
366 b
= (V1
.basis() + V1
.complement().basis())
367 return V
.span_of_basis(b
)
371 def _charpoly_coeff(self
, i
):
373 Return the coefficient polynomial "a_{i}" of this algebra's
374 general characteristic polynomial.
376 Having this be a separate cached method lets us compute and
377 store the trace/determinant (a_{r-1} and a_{0} respectively)
378 separate from the entire characteristic polynomial.
380 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
381 R
= A_of_x
.base_ring()
383 # Guaranteed by theory
386 # Danger: the in-place modification is done for performance
387 # reasons (reconstructing a matrix with huge polynomial
388 # entries is slow), but I don't know how cached_method works,
389 # so it's highly possible that we're modifying some global
390 # list variable by reference, here. In other words, you
391 # probably shouldn't call this method twice on the same
392 # algebra, at the same time, in two threads
393 Ai_orig
= A_of_x
.column(i
)
394 A_of_x
.set_column(i
,xr
)
395 numerator
= A_of_x
.det()
396 A_of_x
.set_column(i
,Ai_orig
)
398 # We're relying on the theory here to ensure that each a_i is
399 # indeed back in R, and the added negative signs are to make
400 # the whole charpoly expression sum to zero.
401 return R(-numerator
/detA
)
405 def _charpoly_matrix_system(self
):
407 Compute the matrix whose entries A_ij are polynomials in
408 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
409 corresponding to `x^r` and the determinent of the matrix A =
410 [A_ij]. In other words, all of the fixed (cachable) data needed
411 to compute the coefficients of the characteristic polynomial.
416 # Construct a new algebra over a multivariate polynomial ring...
417 names
= ['X' + str(i
) for i
in range(1,n
+1)]
418 R
= PolynomialRing(self
.base_ring(), names
)
419 J
= FiniteDimensionalEuclideanJordanAlgebra(R
,
420 self
._multiplication
_table
,
423 idmat
= identity_matrix(J
.base_ring(), n
)
425 W
= self
._charpoly
_basis
_space
()
426 W
= W
.change_ring(R
.fraction_field())
428 # Starting with the standard coordinates x = (X1,X2,...,Xn)
429 # and then converting the entries to W-coordinates allows us
430 # to pass in the standard coordinates to the charpoly and get
431 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
434 # W.coordinates(x^2) eval'd at (standard z-coords)
438 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
440 # We want the middle equivalent thing in our matrix, but use
441 # the first equivalent thing instead so that we can pass in
442 # standard coordinates.
443 x
= J(vector(R
, R
.gens()))
444 l1
= [column_matrix(W
.coordinates((x
**k
).vector())) for k
in range(r
)]
445 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
446 A_of_x
= block_matrix(R
, 1, n
, (l1
+ l2
))
447 xr
= W
.coordinates((x
**r
).vector())
448 return (A_of_x
, x
, xr
, A_of_x
.det())
452 def characteristic_polynomial(self
):
457 This implementation doesn't guarantee that the polynomial
458 denominator in the coefficients is not identically zero, so
459 theoretically it could crash. The way that this is handled
460 in e.g. Faraut and Koranyi is to use a basis that guarantees
461 the denominator is non-zero. But, doing so requires knowledge
462 of at least one regular element, and we don't even know how
463 to do that. The trade-off is that, if we use the standard basis,
464 the resulting polynomial will accept the "usual" coordinates. In
465 other words, we don't have to do a change of basis before e.g.
466 computing the trace or determinant.
470 The characteristic polynomial in the spin algebra is given in
471 Alizadeh, Example 11.11::
473 sage: J = JordanSpinEJA(3)
474 sage: p = J.characteristic_polynomial(); p
475 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
476 sage: xvec = J.one().vector()
484 # The list of coefficient polynomials a_1, a_2, ..., a_n.
485 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
487 # We go to a bit of trouble here to reorder the
488 # indeterminates, so that it's easier to evaluate the
489 # characteristic polynomial at x's coordinates and get back
490 # something in terms of t, which is what we want.
492 S
= PolynomialRing(self
.base_ring(),'t')
494 S
= PolynomialRing(S
, R
.variable_names())
497 # Note: all entries past the rth should be zero. The
498 # coefficient of the highest power (x^r) is 1, but it doesn't
499 # appear in the solution vector which contains coefficients
500 # for the other powers (to make them sum to x^r).
502 a
[r
] = 1 # corresponds to x^r
504 # When the rank is equal to the dimension, trying to
505 # assign a[r] goes out-of-bounds.
506 a
.append(1) # corresponds to x^r
508 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
511 def inner_product(self
, x
, y
):
513 The inner product associated with this Euclidean Jordan algebra.
515 Defaults to the trace inner product, but can be overridden by
516 subclasses if they are sure that the necessary properties are
521 The inner product must satisfy its axiom for this algebra to truly
522 be a Euclidean Jordan Algebra::
524 sage: set_random_seed()
525 sage: J = random_eja()
526 sage: x = J.random_element()
527 sage: y = J.random_element()
528 sage: z = J.random_element()
529 sage: (x*y).inner_product(z) == y.inner_product(x*z)
533 if (not x
in self
) or (not y
in self
):
534 raise TypeError("arguments must live in this algebra")
535 return x
.trace_inner_product(y
)
538 def natural_basis(self
):
540 Return a more-natural representation of this algebra's basis.
542 Every finite-dimensional Euclidean Jordan Algebra is a direct
543 sum of five simple algebras, four of which comprise Hermitian
544 matrices. This method returns the original "natural" basis
545 for our underlying vector space. (Typically, the natural basis
546 is used to construct the multiplication table in the first place.)
548 Note that this will always return a matrix. The standard basis
549 in `R^n` will be returned as `n`-by-`1` column matrices.
553 sage: J = RealSymmetricEJA(2)
556 sage: J.natural_basis()
564 sage: J = JordanSpinEJA(2)
567 sage: J.natural_basis()
574 if self
._natural
_basis
is None:
575 return tuple( b
.vector().column() for b
in self
.basis() )
577 return self
._natural
_basis
582 Return the rank of this EJA.
584 if self
._rank
is None:
585 raise ValueError("no rank specified at genesis")
590 class Element(FiniteDimensionalAlgebraElement
):
592 An element of a Euclidean Jordan algebra.
597 Oh man, I should not be doing this. This hides the "disabled"
598 methods ``left_matrix`` and ``matrix`` from introspection;
599 in particular it removes them from tab-completion.
601 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
602 dir(self
.__class
__) )
605 def __init__(self
, A
, elt
=None):
609 The identity in `S^n` is converted to the identity in the EJA::
611 sage: J = RealSymmetricEJA(3)
612 sage: I = identity_matrix(QQ,3)
613 sage: J(I) == J.one()
616 This skew-symmetric matrix can't be represented in the EJA::
618 sage: J = RealSymmetricEJA(3)
619 sage: A = matrix(QQ,3, lambda i,j: i-j)
621 Traceback (most recent call last):
623 ArithmeticError: vector is not in free module
626 # Goal: if we're given a matrix, and if it lives in our
627 # parent algebra's "natural ambient space," convert it
628 # into an algebra element.
630 # The catch is, we make a recursive call after converting
631 # the given matrix into a vector that lives in the algebra.
632 # This we need to try the parent class initializer first,
633 # to avoid recursing forever if we're given something that
634 # already fits into the algebra, but also happens to live
635 # in the parent's "natural ambient space" (this happens with
638 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, elt
)
640 natural_basis
= A
.natural_basis()
641 if elt
in natural_basis
[0].matrix_space():
642 # Thanks for nothing! Matrix spaces aren't vector
643 # spaces in Sage, so we have to figure out its
644 # natural-basis coordinates ourselves.
645 V
= VectorSpace(elt
.base_ring(), elt
.nrows()**2)
646 W
= V
.span( _mat2vec(s
) for s
in natural_basis
)
647 coords
= W
.coordinates(_mat2vec(elt
))
648 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, coords
)
650 def __pow__(self
, n
):
652 Return ``self`` raised to the power ``n``.
654 Jordan algebras are always power-associative; see for
655 example Faraut and Koranyi, Proposition II.1.2 (ii).
659 We have to override this because our superclass uses row vectors
660 instead of column vectors! We, on the other hand, assume column
665 sage: set_random_seed()
666 sage: x = random_eja().random_element()
667 sage: x.operator_matrix()*x.vector() == (x^2).vector()
670 A few examples of power-associativity::
672 sage: set_random_seed()
673 sage: x = random_eja().random_element()
674 sage: x*(x*x)*(x*x) == x^5
676 sage: (x*x)*(x*x*x) == x^5
679 We also know that powers operator-commute (Koecher, Chapter
682 sage: set_random_seed()
683 sage: x = random_eja().random_element()
684 sage: m = ZZ.random_element(0,10)
685 sage: n = ZZ.random_element(0,10)
686 sage: Lxm = (x^m).operator_matrix()
687 sage: Lxn = (x^n).operator_matrix()
688 sage: Lxm*Lxn == Lxn*Lxm
698 return A( (self
.operator_matrix()**(n
-1))*self
.vector() )
701 def apply_univariate_polynomial(self
, p
):
703 Apply the univariate polynomial ``p`` to this element.
705 A priori, SageMath won't allow us to apply a univariate
706 polynomial to an element of an EJA, because we don't know
707 that EJAs are rings (they are usually not associative). Of
708 course, we know that EJAs are power-associative, so the
709 operation is ultimately kosher. This function sidesteps
710 the CAS to get the answer we want and expect.
714 sage: R = PolynomialRing(QQ, 't')
716 sage: p = t^4 - t^3 + 5*t - 2
717 sage: J = RealCartesianProductEJA(5)
718 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
723 We should always get back an element of the algebra::
725 sage: set_random_seed()
726 sage: p = PolynomialRing(QQ, 't').random_element()
727 sage: J = random_eja()
728 sage: x = J.random_element()
729 sage: x.apply_univariate_polynomial(p) in J
733 if len(p
.variables()) > 1:
734 raise ValueError("not a univariate polynomial")
737 # Convert the coeficcients to the parent's base ring,
738 # because a priori they might live in an (unnecessarily)
739 # larger ring for which P.sum() would fail below.
740 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
741 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
744 def characteristic_polynomial(self
):
746 Return the characteristic polynomial of this element.
750 The rank of `R^3` is three, and the minimal polynomial of
751 the identity element is `(t-1)` from which it follows that
752 the characteristic polynomial should be `(t-1)^3`::
754 sage: J = RealCartesianProductEJA(3)
755 sage: J.one().characteristic_polynomial()
756 t^3 - 3*t^2 + 3*t - 1
758 Likewise, the characteristic of the zero element in the
759 rank-three algebra `R^{n}` should be `t^{3}`::
761 sage: J = RealCartesianProductEJA(3)
762 sage: J.zero().characteristic_polynomial()
765 The characteristic polynomial of an element should evaluate
766 to zero on that element::
768 sage: set_random_seed()
769 sage: x = RealCartesianProductEJA(3).random_element()
770 sage: p = x.characteristic_polynomial()
771 sage: x.apply_univariate_polynomial(p)
775 p
= self
.parent().characteristic_polynomial()
776 return p(*self
.vector())
779 def inner_product(self
, other
):
781 Return the parent algebra's inner product of myself and ``other``.
785 The inner product in the Jordan spin algebra is the usual
786 inner product on `R^n` (this example only works because the
787 basis for the Jordan algebra is the standard basis in `R^n`)::
789 sage: J = JordanSpinEJA(3)
790 sage: x = vector(QQ,[1,2,3])
791 sage: y = vector(QQ,[4,5,6])
792 sage: x.inner_product(y)
794 sage: J(x).inner_product(J(y))
797 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
798 multiplication is the usual matrix multiplication in `S^n`,
799 so the inner product of the identity matrix with itself
802 sage: J = RealSymmetricEJA(3)
803 sage: J.one().inner_product(J.one())
806 Likewise, the inner product on `C^n` is `<X,Y> =
807 Re(trace(X*Y))`, where we must necessarily take the real
808 part because the product of Hermitian matrices may not be
811 sage: J = ComplexHermitianEJA(3)
812 sage: J.one().inner_product(J.one())
815 Ditto for the quaternions::
817 sage: J = QuaternionHermitianEJA(3)
818 sage: J.one().inner_product(J.one())
823 Ensure that we can always compute an inner product, and that
824 it gives us back a real number::
826 sage: set_random_seed()
827 sage: J = random_eja()
828 sage: x = J.random_element()
829 sage: y = J.random_element()
830 sage: x.inner_product(y) in RR
836 raise TypeError("'other' must live in the same algebra")
838 return P
.inner_product(self
, other
)
841 def operator_commutes_with(self
, other
):
843 Return whether or not this element operator-commutes
848 The definition of a Jordan algebra says that any element
849 operator-commutes with its square::
851 sage: set_random_seed()
852 sage: x = random_eja().random_element()
853 sage: x.operator_commutes_with(x^2)
858 Test Lemma 1 from Chapter III of Koecher::
860 sage: set_random_seed()
861 sage: J = random_eja()
862 sage: u = J.random_element()
863 sage: v = J.random_element()
864 sage: lhs = u.operator_commutes_with(u*v)
865 sage: rhs = v.operator_commutes_with(u^2)
870 if not other
in self
.parent():
871 raise TypeError("'other' must live in the same algebra")
873 A
= self
.operator_matrix()
874 B
= other
.operator_matrix()
880 Return my determinant, the product of my eigenvalues.
884 sage: J = JordanSpinEJA(2)
885 sage: e0,e1 = J.gens()
886 sage: x = sum( J.gens() )
892 sage: J = JordanSpinEJA(3)
893 sage: e0,e1,e2 = J.gens()
894 sage: x = sum( J.gens() )
900 An element is invertible if and only if its determinant is
903 sage: set_random_seed()
904 sage: x = random_eja().random_element()
905 sage: x.is_invertible() == (x.det() != 0)
911 p
= P
._charpoly
_coeff
(0)
912 # The _charpoly_coeff function already adds the factor of
913 # -1 to ensure that _charpoly_coeff(0) is really what
914 # appears in front of t^{0} in the charpoly. However,
915 # we want (-1)^r times THAT for the determinant.
916 return ((-1)**r
)*p(*self
.vector())
921 Return the Jordan-multiplicative inverse of this element.
925 We appeal to the quadratic representation as in Koecher's
926 Theorem 12 in Chapter III, Section 5.
930 The inverse in the spin factor algebra is given in Alizadeh's
933 sage: set_random_seed()
934 sage: n = ZZ.random_element(1,10)
935 sage: J = JordanSpinEJA(n)
936 sage: x = J.random_element()
937 sage: while x.is_zero():
938 ....: x = J.random_element()
939 sage: x_vec = x.vector()
941 sage: x_bar = x_vec[1:]
942 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
943 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
944 sage: x_inverse = coeff*inv_vec
945 sage: x.inverse() == J(x_inverse)
950 The identity element is its own inverse::
952 sage: set_random_seed()
953 sage: J = random_eja()
954 sage: J.one().inverse() == J.one()
957 If an element has an inverse, it acts like one::
959 sage: set_random_seed()
960 sage: J = random_eja()
961 sage: x = J.random_element()
962 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
965 The inverse of the inverse is what we started with::
967 sage: set_random_seed()
968 sage: J = random_eja()
969 sage: x = J.random_element()
970 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
973 The zero element is never invertible::
975 sage: set_random_seed()
976 sage: J = random_eja().zero().inverse()
977 Traceback (most recent call last):
979 ValueError: element is not invertible
982 if not self
.is_invertible():
983 raise ValueError("element is not invertible")
986 return P(self
.quadratic_representation().inverse()*self
.vector())
989 def is_invertible(self
):
991 Return whether or not this element is invertible.
993 We can't use the superclass method because it relies on
994 the algebra being associative.
998 The usual way to do this is to check if the determinant is
999 zero, but we need the characteristic polynomial for the
1000 determinant. The minimal polynomial is a lot easier to get,
1001 so we use Corollary 2 in Chapter V of Koecher to check
1002 whether or not the paren't algebra's zero element is a root
1003 of this element's minimal polynomial.
1007 The identity element is always invertible::
1009 sage: set_random_seed()
1010 sage: J = random_eja()
1011 sage: J.one().is_invertible()
1014 The zero element is never invertible::
1016 sage: set_random_seed()
1017 sage: J = random_eja()
1018 sage: J.zero().is_invertible()
1022 zero
= self
.parent().zero()
1023 p
= self
.minimal_polynomial()
1024 return not (p(zero
) == zero
)
1027 def is_nilpotent(self
):
1029 Return whether or not some power of this element is zero.
1031 The superclass method won't work unless we're in an
1032 associative algebra, and we aren't. However, we generate
1033 an assocoative subalgebra and we're nilpotent there if and
1034 only if we're nilpotent here (probably).
1038 The identity element is never nilpotent::
1040 sage: set_random_seed()
1041 sage: random_eja().one().is_nilpotent()
1044 The additive identity is always nilpotent::
1046 sage: set_random_seed()
1047 sage: random_eja().zero().is_nilpotent()
1051 # The element we're going to call "is_nilpotent()" on.
1052 # Either myself, interpreted as an element of a finite-
1053 # dimensional algebra, or an element of an associative
1057 if self
.parent().is_associative():
1058 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1060 V
= self
.span_of_powers()
1061 assoc_subalg
= self
.subalgebra_generated_by()
1062 # Mis-design warning: the basis used for span_of_powers()
1063 # and subalgebra_generated_by() must be the same, and in
1065 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1067 # Recursive call, but should work since elt lives in an
1068 # associative algebra.
1069 return elt
.is_nilpotent()
1072 def is_regular(self
):
1074 Return whether or not this is a regular element.
1078 The identity element always has degree one, but any element
1079 linearly-independent from it is regular::
1081 sage: J = JordanSpinEJA(5)
1082 sage: J.one().is_regular()
1084 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
1085 sage: for x in J.gens():
1086 ....: (J.one() + x).is_regular()
1094 return self
.degree() == self
.parent().rank()
1099 Compute the degree of this element the straightforward way
1100 according to the definition; by appending powers to a list
1101 and figuring out its dimension (that is, whether or not
1102 they're linearly dependent).
1106 sage: J = JordanSpinEJA(4)
1107 sage: J.one().degree()
1109 sage: e0,e1,e2,e3 = J.gens()
1110 sage: (e0 - e1).degree()
1113 In the spin factor algebra (of rank two), all elements that
1114 aren't multiples of the identity are regular::
1116 sage: set_random_seed()
1117 sage: n = ZZ.random_element(1,10)
1118 sage: J = JordanSpinEJA(n)
1119 sage: x = J.random_element()
1120 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
1124 return self
.span_of_powers().dimension()
1127 def left_matrix(self
):
1129 Our parent class defines ``left_matrix`` and ``matrix``
1130 methods whose names are misleading. We don't want them.
1132 raise NotImplementedError("use operator_matrix() instead")
1134 matrix
= left_matrix
1137 def minimal_polynomial(self
):
1139 Return the minimal polynomial of this element,
1140 as a function of the variable `t`.
1144 We restrict ourselves to the associative subalgebra
1145 generated by this element, and then return the minimal
1146 polynomial of this element's operator matrix (in that
1147 subalgebra). This works by Baes Proposition 2.3.16.
1151 The minimal polynomial of the identity and zero elements are
1154 sage: set_random_seed()
1155 sage: J = random_eja()
1156 sage: J.one().minimal_polynomial()
1158 sage: J.zero().minimal_polynomial()
1161 The degree of an element is (by one definition) the degree
1162 of its minimal polynomial::
1164 sage: set_random_seed()
1165 sage: x = random_eja().random_element()
1166 sage: x.degree() == x.minimal_polynomial().degree()
1169 The minimal polynomial and the characteristic polynomial coincide
1170 and are known (see Alizadeh, Example 11.11) for all elements of
1171 the spin factor algebra that aren't scalar multiples of the
1174 sage: set_random_seed()
1175 sage: n = ZZ.random_element(2,10)
1176 sage: J = JordanSpinEJA(n)
1177 sage: y = J.random_element()
1178 sage: while y == y.coefficient(0)*J.one():
1179 ....: y = J.random_element()
1180 sage: y0 = y.vector()[0]
1181 sage: y_bar = y.vector()[1:]
1182 sage: actual = y.minimal_polynomial()
1183 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1184 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1185 sage: bool(actual == expected)
1188 The minimal polynomial should always kill its element::
1190 sage: set_random_seed()
1191 sage: x = random_eja().random_element()
1192 sage: p = x.minimal_polynomial()
1193 sage: x.apply_univariate_polynomial(p)
1197 V
= self
.span_of_powers()
1198 assoc_subalg
= self
.subalgebra_generated_by()
1199 # Mis-design warning: the basis used for span_of_powers()
1200 # and subalgebra_generated_by() must be the same, and in
1202 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1204 # We get back a symbolic polynomial in 'x' but want a real
1205 # polynomial in 't'.
1206 p_of_x
= elt
.operator_matrix().minimal_polynomial()
1207 return p_of_x
.change_variable_name('t')
1210 def natural_representation(self
):
1212 Return a more-natural representation of this element.
1214 Every finite-dimensional Euclidean Jordan Algebra is a
1215 direct sum of five simple algebras, four of which comprise
1216 Hermitian matrices. This method returns the original
1217 "natural" representation of this element as a Hermitian
1218 matrix, if it has one. If not, you get the usual representation.
1222 sage: J = ComplexHermitianEJA(3)
1225 sage: J.one().natural_representation()
1235 sage: J = QuaternionHermitianEJA(3)
1238 sage: J.one().natural_representation()
1239 [1 0 0 0 0 0 0 0 0 0 0 0]
1240 [0 1 0 0 0 0 0 0 0 0 0 0]
1241 [0 0 1 0 0 0 0 0 0 0 0 0]
1242 [0 0 0 1 0 0 0 0 0 0 0 0]
1243 [0 0 0 0 1 0 0 0 0 0 0 0]
1244 [0 0 0 0 0 1 0 0 0 0 0 0]
1245 [0 0 0 0 0 0 1 0 0 0 0 0]
1246 [0 0 0 0 0 0 0 1 0 0 0 0]
1247 [0 0 0 0 0 0 0 0 1 0 0 0]
1248 [0 0 0 0 0 0 0 0 0 1 0 0]
1249 [0 0 0 0 0 0 0 0 0 0 1 0]
1250 [0 0 0 0 0 0 0 0 0 0 0 1]
1253 B
= self
.parent().natural_basis()
1254 W
= B
[0].matrix_space()
1255 return W
.linear_combination(zip(self
.vector(), B
))
1260 Return the left-multiplication-by-this-element
1261 operator on the ambient algebra.
1265 sage: set_random_seed()
1266 sage: J = random_eja()
1267 sage: x = J.random_element()
1268 sage: y = J.random_element()
1269 sage: x.operator()(y) == x*y
1271 sage: y.operator()(x) == x*y
1276 return FiniteDimensionalEuclideanJordanAlgebraMorphism(
1278 self
.operator_matrix() )
1282 def operator_matrix(self
):
1284 Return the matrix that represents left- (or right-)
1285 multiplication by this element in the parent algebra.
1287 We implement this ourselves to work around the fact that
1288 our parent class represents everything with row vectors.
1292 Test the first polarization identity from my notes, Koecher Chapter
1293 III, or from Baes (2.3)::
1295 sage: set_random_seed()
1296 sage: J = random_eja()
1297 sage: x = J.random_element()
1298 sage: y = J.random_element()
1299 sage: Lx = x.operator_matrix()
1300 sage: Ly = y.operator_matrix()
1301 sage: Lxx = (x*x).operator_matrix()
1302 sage: Lxy = (x*y).operator_matrix()
1303 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
1306 Test the second polarization identity from my notes or from
1309 sage: set_random_seed()
1310 sage: J = random_eja()
1311 sage: x = J.random_element()
1312 sage: y = J.random_element()
1313 sage: z = J.random_element()
1314 sage: Lx = x.operator_matrix()
1315 sage: Ly = y.operator_matrix()
1316 sage: Lz = z.operator_matrix()
1317 sage: Lzy = (z*y).operator_matrix()
1318 sage: Lxy = (x*y).operator_matrix()
1319 sage: Lxz = (x*z).operator_matrix()
1320 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
1323 Test the third polarization identity from my notes or from
1326 sage: set_random_seed()
1327 sage: J = random_eja()
1328 sage: u = J.random_element()
1329 sage: y = J.random_element()
1330 sage: z = J.random_element()
1331 sage: Lu = u.operator_matrix()
1332 sage: Ly = y.operator_matrix()
1333 sage: Lz = z.operator_matrix()
1334 sage: Lzy = (z*y).operator_matrix()
1335 sage: Luy = (u*y).operator_matrix()
1336 sage: Luz = (u*z).operator_matrix()
1337 sage: Luyz = (u*(y*z)).operator_matrix()
1338 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
1339 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
1340 sage: bool(lhs == rhs)
1344 fda_elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
1345 return fda_elt
.matrix().transpose()
1348 def quadratic_representation(self
, other
=None):
1350 Return the quadratic representation of this element.
1354 The explicit form in the spin factor algebra is given by
1355 Alizadeh's Example 11.12::
1357 sage: set_random_seed()
1358 sage: n = ZZ.random_element(1,10)
1359 sage: J = JordanSpinEJA(n)
1360 sage: x = J.random_element()
1361 sage: x_vec = x.vector()
1363 sage: x_bar = x_vec[1:]
1364 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1365 sage: B = 2*x0*x_bar.row()
1366 sage: C = 2*x0*x_bar.column()
1367 sage: D = identity_matrix(QQ, n-1)
1368 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1369 sage: D = D + 2*x_bar.tensor_product(x_bar)
1370 sage: Q = block_matrix(2,2,[A,B,C,D])
1371 sage: Q == x.quadratic_representation().operator_matrix()
1374 Test all of the properties from Theorem 11.2 in Alizadeh::
1376 sage: set_random_seed()
1377 sage: J = random_eja()
1378 sage: x = J.random_element()
1379 sage: y = J.random_element()
1380 sage: Lx = x.operator()
1381 sage: Lxx = (x*x).operator()
1382 sage: Qx = x.quadratic_representation()
1383 sage: Qy = y.quadratic_representation()
1384 sage: Qxy = x.quadratic_representation(y)
1385 sage: Qex = J.one().quadratic_representation(x)
1386 sage: n = ZZ.random_element(10)
1387 sage: Qxn = (x^n).quadratic_representation()
1391 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1396 sage: alpha = QQ.random_element()
1397 sage: (alpha*x).quadratic_representation() == (alpha^2)*Qx
1402 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1405 sage: not x.is_invertible() or (
1408 ....: x.inverse().quadratic_representation() )
1411 sage: Qxy(J.one()) == x*y
1416 sage: not x.is_invertible() or (
1417 ....: x.quadratic_representation(x.inverse())*Qx
1418 ....: == Qx*x.quadratic_representation(x.inverse()) )
1421 sage: not x.is_invertible() or (
1422 ....: x.quadratic_representation(x.inverse())*Qx
1424 ....: 2*x.operator()*Qex - Qx )
1427 sage: 2*x.operator()*Qex - Qx == Lxx
1432 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1442 sage: not x.is_invertible() or (
1443 ....: Qx*x.inverse().operator() == Lx )
1448 sage: not x.operator_commutes_with(y) or (
1449 ....: Qx(y)^n == Qxn(y^n) )
1455 elif not other
in self
.parent():
1456 raise TypeError("'other' must live in the same algebra")
1459 M
= other
.operator()
1460 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1463 def span_of_powers(self
):
1465 Return the vector space spanned by successive powers of
1468 # The dimension of the subalgebra can't be greater than
1469 # the big algebra, so just put everything into a list
1470 # and let span() get rid of the excess.
1472 # We do the extra ambient_vector_space() in case we're messing
1473 # with polynomials and the direct parent is a module.
1474 V
= self
.vector().parent().ambient_vector_space()
1475 return V
.span( (self
**d
).vector() for d
in xrange(V
.dimension()) )
1478 def subalgebra_generated_by(self
):
1480 Return the associative subalgebra of the parent EJA generated
1485 sage: set_random_seed()
1486 sage: x = random_eja().random_element()
1487 sage: x.subalgebra_generated_by().is_associative()
1490 Squaring in the subalgebra should be the same thing as
1491 squaring in the superalgebra::
1493 sage: set_random_seed()
1494 sage: x = random_eja().random_element()
1495 sage: u = x.subalgebra_generated_by().random_element()
1496 sage: u.operator_matrix()*u.vector() == (u**2).vector()
1500 # First get the subspace spanned by the powers of myself...
1501 V
= self
.span_of_powers()
1502 F
= self
.base_ring()
1504 # Now figure out the entries of the right-multiplication
1505 # matrix for the successive basis elements b0, b1,... of
1508 for b_right
in V
.basis():
1509 eja_b_right
= self
.parent()(b_right
)
1511 # The first row of the right-multiplication matrix by
1512 # b1 is what we get if we apply that matrix to b1. The
1513 # second row of the right multiplication matrix by b1
1514 # is what we get when we apply that matrix to b2...
1516 # IMPORTANT: this assumes that all vectors are COLUMN
1517 # vectors, unlike our superclass (which uses row vectors).
1518 for b_left
in V
.basis():
1519 eja_b_left
= self
.parent()(b_left
)
1520 # Multiply in the original EJA, but then get the
1521 # coordinates from the subalgebra in terms of its
1523 this_row
= V
.coordinates((eja_b_left
*eja_b_right
).vector())
1524 b_right_rows
.append(this_row
)
1525 b_right_matrix
= matrix(F
, b_right_rows
)
1526 mats
.append(b_right_matrix
)
1528 # It's an algebra of polynomials in one element, and EJAs
1529 # are power-associative.
1531 # TODO: choose generator names intelligently.
1532 return FiniteDimensionalEuclideanJordanAlgebra(F
, mats
, assume_associative
=True, names
='f')
1535 def subalgebra_idempotent(self
):
1537 Find an idempotent in the associative subalgebra I generate
1538 using Proposition 2.3.5 in Baes.
1542 sage: set_random_seed()
1543 sage: J = random_eja()
1544 sage: x = J.random_element()
1545 sage: while x.is_nilpotent():
1546 ....: x = J.random_element()
1547 sage: c = x.subalgebra_idempotent()
1552 if self
.is_nilpotent():
1553 raise ValueError("this only works with non-nilpotent elements!")
1555 V
= self
.span_of_powers()
1556 J
= self
.subalgebra_generated_by()
1557 # Mis-design warning: the basis used for span_of_powers()
1558 # and subalgebra_generated_by() must be the same, and in
1560 u
= J(V
.coordinates(self
.vector()))
1562 # The image of the matrix of left-u^m-multiplication
1563 # will be minimal for some natural number s...
1565 minimal_dim
= V
.dimension()
1566 for i
in xrange(1, V
.dimension()):
1567 this_dim
= (u
**i
).operator_matrix().image().dimension()
1568 if this_dim
< minimal_dim
:
1569 minimal_dim
= this_dim
1572 # Now minimal_matrix should correspond to the smallest
1573 # non-zero subspace in Baes's (or really, Koecher's)
1576 # However, we need to restrict the matrix to work on the
1577 # subspace... or do we? Can't we just solve, knowing that
1578 # A(c) = u^(s+1) should have a solution in the big space,
1581 # Beware, solve_right() means that we're using COLUMN vectors.
1582 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1584 A
= u_next
.operator_matrix()
1585 c_coordinates
= A
.solve_right(u_next
.vector())
1587 # Now c_coordinates is the idempotent we want, but it's in
1588 # the coordinate system of the subalgebra.
1590 # We need the basis for J, but as elements of the parent algebra.
1592 basis
= [self
.parent(v
) for v
in V
.basis()]
1593 return self
.parent().linear_combination(zip(c_coordinates
, basis
))
1598 Return my trace, the sum of my eigenvalues.
1602 sage: J = JordanSpinEJA(3)
1603 sage: x = sum(J.gens())
1609 sage: J = RealCartesianProductEJA(5)
1610 sage: J.one().trace()
1615 The trace of an element is a real number::
1617 sage: set_random_seed()
1618 sage: J = random_eja()
1619 sage: J.random_element().trace() in J.base_ring()
1625 p
= P
._charpoly
_coeff
(r
-1)
1626 # The _charpoly_coeff function already adds the factor of
1627 # -1 to ensure that _charpoly_coeff(r-1) is really what
1628 # appears in front of t^{r-1} in the charpoly. However,
1629 # we want the negative of THAT for the trace.
1630 return -p(*self
.vector())
1633 def trace_inner_product(self
, other
):
1635 Return the trace inner product of myself and ``other``.
1639 The trace inner product is commutative::
1641 sage: set_random_seed()
1642 sage: J = random_eja()
1643 sage: x = J.random_element(); y = J.random_element()
1644 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1647 The trace inner product is bilinear::
1649 sage: set_random_seed()
1650 sage: J = random_eja()
1651 sage: x = J.random_element()
1652 sage: y = J.random_element()
1653 sage: z = J.random_element()
1654 sage: a = QQ.random_element();
1655 sage: actual = (a*(x+z)).trace_inner_product(y)
1656 sage: expected = ( a*x.trace_inner_product(y) +
1657 ....: a*z.trace_inner_product(y) )
1658 sage: actual == expected
1660 sage: actual = x.trace_inner_product(a*(y+z))
1661 sage: expected = ( a*x.trace_inner_product(y) +
1662 ....: a*x.trace_inner_product(z) )
1663 sage: actual == expected
1666 The trace inner product satisfies the compatibility
1667 condition in the definition of a Euclidean Jordan algebra::
1669 sage: set_random_seed()
1670 sage: J = random_eja()
1671 sage: x = J.random_element()
1672 sage: y = J.random_element()
1673 sage: z = J.random_element()
1674 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1678 if not other
in self
.parent():
1679 raise TypeError("'other' must live in the same algebra")
1681 return (self
*other
).trace()
1684 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1686 Return the Euclidean Jordan Algebra corresponding to the set
1687 `R^n` under the Hadamard product.
1689 Note: this is nothing more than the Cartesian product of ``n``
1690 copies of the spin algebra. Once Cartesian product algebras
1691 are implemented, this can go.
1695 This multiplication table can be verified by hand::
1697 sage: J = RealCartesianProductEJA(3)
1698 sage: e0,e1,e2 = J.gens()
1714 def __classcall_private__(cls
, n
, field
=QQ
):
1715 # The FiniteDimensionalAlgebra constructor takes a list of
1716 # matrices, the ith representing right multiplication by the ith
1717 # basis element in the vector space. So if e_1 = (1,0,0), then
1718 # right (Hadamard) multiplication of x by e_1 picks out the first
1719 # component of x; and likewise for the ith basis element e_i.
1720 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
1721 for i
in xrange(n
) ]
1723 fdeja
= super(RealCartesianProductEJA
, cls
)
1724 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=n
)
1726 def inner_product(self
, x
, y
):
1727 return _usual_ip(x
,y
)
1732 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1736 For now, we choose a random natural number ``n`` (greater than zero)
1737 and then give you back one of the following:
1739 * The cartesian product of the rational numbers ``n`` times; this is
1740 ``QQ^n`` with the Hadamard product.
1742 * The Jordan spin algebra on ``QQ^n``.
1744 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1747 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1748 in the space of ``2n``-by-``2n`` real symmetric matrices.
1750 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1751 in the space of ``4n``-by-``4n`` real symmetric matrices.
1753 Later this might be extended to return Cartesian products of the
1759 Euclidean Jordan algebra of degree...
1763 # The max_n component lets us choose different upper bounds on the
1764 # value "n" that gets passed to the constructor. This is needed
1765 # because e.g. R^{10} is reasonable to test, while the Hermitian
1766 # 10-by-10 quaternion matrices are not.
1767 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
1769 (RealSymmetricEJA
, 5),
1770 (ComplexHermitianEJA
, 4),
1771 (QuaternionHermitianEJA
, 3)])
1772 n
= ZZ
.random_element(1, max_n
)
1773 return constructor(n
, field
=QQ
)
1777 def _real_symmetric_basis(n
, field
=QQ
):
1779 Return a basis for the space of real symmetric n-by-n matrices.
1781 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1785 for j
in xrange(i
+1):
1786 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1790 # Beware, orthogonal but not normalized!
1791 Sij
= Eij
+ Eij
.transpose()
1796 def _complex_hermitian_basis(n
, field
=QQ
):
1798 Returns a basis for the space of complex Hermitian n-by-n matrices.
1802 sage: set_random_seed()
1803 sage: n = ZZ.random_element(1,5)
1804 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1808 F
= QuadraticField(-1, 'I')
1811 # This is like the symmetric case, but we need to be careful:
1813 # * We want conjugate-symmetry, not just symmetry.
1814 # * The diagonal will (as a result) be real.
1818 for j
in xrange(i
+1):
1819 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1821 Sij
= _embed_complex_matrix(Eij
)
1824 # Beware, orthogonal but not normalized! The second one
1825 # has a minus because it's conjugated.
1826 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
1828 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
1833 def _quaternion_hermitian_basis(n
, field
=QQ
):
1835 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1839 sage: set_random_seed()
1840 sage: n = ZZ.random_element(1,5)
1841 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1845 Q
= QuaternionAlgebra(QQ
,-1,-1)
1848 # This is like the symmetric case, but we need to be careful:
1850 # * We want conjugate-symmetry, not just symmetry.
1851 # * The diagonal will (as a result) be real.
1855 for j
in xrange(i
+1):
1856 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1858 Sij
= _embed_quaternion_matrix(Eij
)
1861 # Beware, orthogonal but not normalized! The second,
1862 # third, and fourth ones have a minus because they're
1864 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
1866 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
1868 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
1870 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
1876 return vector(m
.base_ring(), m
.list())
1879 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
1881 def _multiplication_table_from_matrix_basis(basis
):
1883 At least three of the five simple Euclidean Jordan algebras have the
1884 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1885 multiplication on the right is matrix multiplication. Given a basis
1886 for the underlying matrix space, this function returns a
1887 multiplication table (obtained by looping through the basis
1888 elements) for an algebra of those matrices. A reordered copy
1889 of the basis is also returned to work around the fact that
1890 the ``span()`` in this function will change the order of the basis
1891 from what we think it is, to... something else.
1893 # In S^2, for example, we nominally have four coordinates even
1894 # though the space is of dimension three only. The vector space V
1895 # is supposed to hold the entire long vector, and the subspace W
1896 # of V will be spanned by the vectors that arise from symmetric
1897 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1898 field
= basis
[0].base_ring()
1899 dimension
= basis
[0].nrows()
1901 V
= VectorSpace(field
, dimension
**2)
1902 W
= V
.span( _mat2vec(s
) for s
in basis
)
1904 # Taking the span above reorders our basis (thanks, jerk!) so we
1905 # need to put our "matrix basis" in the same order as the
1906 # (reordered) vector basis.
1907 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
1911 # Brute force the multiplication-by-s matrix by looping
1912 # through all elements of the basis and doing the computation
1913 # to find out what the corresponding row should be. BEWARE:
1914 # these multiplication tables won't be symmetric! It therefore
1915 # becomes REALLY IMPORTANT that the underlying algebra
1916 # constructor uses ROW vectors and not COLUMN vectors. That's
1917 # why we're computing rows here and not columns.
1920 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
1921 Q_rows
.append(W
.coordinates(this_row
))
1922 Q
= matrix(field
, W
.dimension(), Q_rows
)
1928 def _embed_complex_matrix(M
):
1930 Embed the n-by-n complex matrix ``M`` into the space of real
1931 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1932 bi` to the block matrix ``[[a,b],[-b,a]]``.
1936 sage: F = QuadraticField(-1,'i')
1937 sage: x1 = F(4 - 2*i)
1938 sage: x2 = F(1 + 2*i)
1941 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1942 sage: _embed_complex_matrix(M)
1951 Embedding is a homomorphism (isomorphism, in fact)::
1953 sage: set_random_seed()
1954 sage: n = ZZ.random_element(5)
1955 sage: F = QuadraticField(-1, 'i')
1956 sage: X = random_matrix(F, n)
1957 sage: Y = random_matrix(F, n)
1958 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1959 sage: expected = _embed_complex_matrix(X*Y)
1960 sage: actual == expected
1966 raise ValueError("the matrix 'M' must be square")
1967 field
= M
.base_ring()
1972 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1974 # We can drop the imaginaries here.
1975 return block_matrix(field
.base_ring(), n
, blocks
)
1978 def _unembed_complex_matrix(M
):
1980 The inverse of _embed_complex_matrix().
1984 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1985 ....: [-2, 1, -4, 3],
1986 ....: [ 9, 10, 11, 12],
1987 ....: [-10, 9, -12, 11] ])
1988 sage: _unembed_complex_matrix(A)
1990 [ 10*i + 9 12*i + 11]
1994 Unembedding is the inverse of embedding::
1996 sage: set_random_seed()
1997 sage: F = QuadraticField(-1, 'i')
1998 sage: M = random_matrix(F, 3)
1999 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
2005 raise ValueError("the matrix 'M' must be square")
2006 if not n
.mod(2).is_zero():
2007 raise ValueError("the matrix 'M' must be a complex embedding")
2009 F
= QuadraticField(-1, 'i')
2012 # Go top-left to bottom-right (reading order), converting every
2013 # 2-by-2 block we see to a single complex element.
2015 for k
in xrange(n
/2):
2016 for j
in xrange(n
/2):
2017 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
2018 if submat
[0,0] != submat
[1,1]:
2019 raise ValueError('bad on-diagonal submatrix')
2020 if submat
[0,1] != -submat
[1,0]:
2021 raise ValueError('bad off-diagonal submatrix')
2022 z
= submat
[0,0] + submat
[0,1]*i
2025 return matrix(F
, n
/2, elements
)
2028 def _embed_quaternion_matrix(M
):
2030 Embed the n-by-n quaternion matrix ``M`` into the space of real
2031 matrices of size 4n-by-4n by first sending each quaternion entry
2032 `z = a + bi + cj + dk` to the block-complex matrix
2033 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
2038 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2039 sage: i,j,k = Q.gens()
2040 sage: x = 1 + 2*i + 3*j + 4*k
2041 sage: M = matrix(Q, 1, [[x]])
2042 sage: _embed_quaternion_matrix(M)
2048 Embedding is a homomorphism (isomorphism, in fact)::
2050 sage: set_random_seed()
2051 sage: n = ZZ.random_element(5)
2052 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2053 sage: X = random_matrix(Q, n)
2054 sage: Y = random_matrix(Q, n)
2055 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
2056 sage: expected = _embed_quaternion_matrix(X*Y)
2057 sage: actual == expected
2061 quaternions
= M
.base_ring()
2064 raise ValueError("the matrix 'M' must be square")
2066 F
= QuadraticField(-1, 'i')
2071 t
= z
.coefficient_tuple()
2076 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2077 [-c
+ d
*i
, a
- b
*i
]])
2078 blocks
.append(_embed_complex_matrix(cplx_matrix
))
2080 # We should have real entries by now, so use the realest field
2081 # we've got for the return value.
2082 return block_matrix(quaternions
.base_ring(), n
, blocks
)
2085 def _unembed_quaternion_matrix(M
):
2087 The inverse of _embed_quaternion_matrix().
2091 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2092 ....: [-2, 1, -4, 3],
2093 ....: [-3, 4, 1, -2],
2094 ....: [-4, -3, 2, 1]])
2095 sage: _unembed_quaternion_matrix(M)
2096 [1 + 2*i + 3*j + 4*k]
2100 Unembedding is the inverse of embedding::
2102 sage: set_random_seed()
2103 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2104 sage: M = random_matrix(Q, 3)
2105 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
2111 raise ValueError("the matrix 'M' must be square")
2112 if not n
.mod(4).is_zero():
2113 raise ValueError("the matrix 'M' must be a complex embedding")
2115 Q
= QuaternionAlgebra(QQ
,-1,-1)
2118 # Go top-left to bottom-right (reading order), converting every
2119 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2122 for l
in xrange(n
/4):
2123 for m
in xrange(n
/4):
2124 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
2125 if submat
[0,0] != submat
[1,1].conjugate():
2126 raise ValueError('bad on-diagonal submatrix')
2127 if submat
[0,1] != -submat
[1,0].conjugate():
2128 raise ValueError('bad off-diagonal submatrix')
2129 z
= submat
[0,0].real() + submat
[0,0].imag()*i
2130 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
2133 return matrix(Q
, n
/4, elements
)
2136 # The usual inner product on R^n.
2138 return x
.vector().inner_product(y
.vector())
2140 # The inner product used for the real symmetric simple EJA.
2141 # We keep it as a separate function because e.g. the complex
2142 # algebra uses the same inner product, except divided by 2.
2143 def _matrix_ip(X
,Y
):
2144 X_mat
= X
.natural_representation()
2145 Y_mat
= Y
.natural_representation()
2146 return (X_mat
*Y_mat
).trace()
2149 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2151 The rank-n simple EJA consisting of real symmetric n-by-n
2152 matrices, the usual symmetric Jordan product, and the trace inner
2153 product. It has dimension `(n^2 + n)/2` over the reals.
2157 sage: J = RealSymmetricEJA(2)
2158 sage: e0, e1, e2 = J.gens()
2168 The degree of this algebra is `(n^2 + n) / 2`::
2170 sage: set_random_seed()
2171 sage: n = ZZ.random_element(1,5)
2172 sage: J = RealSymmetricEJA(n)
2173 sage: J.degree() == (n^2 + n)/2
2176 The Jordan multiplication is what we think it is::
2178 sage: set_random_seed()
2179 sage: n = ZZ.random_element(1,5)
2180 sage: J = RealSymmetricEJA(n)
2181 sage: x = J.random_element()
2182 sage: y = J.random_element()
2183 sage: actual = (x*y).natural_representation()
2184 sage: X = x.natural_representation()
2185 sage: Y = y.natural_representation()
2186 sage: expected = (X*Y + Y*X)/2
2187 sage: actual == expected
2189 sage: J(expected) == x*y
2194 def __classcall_private__(cls
, n
, field
=QQ
):
2195 S
= _real_symmetric_basis(n
, field
=field
)
2196 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2198 fdeja
= super(RealSymmetricEJA
, cls
)
2199 return fdeja
.__classcall
_private
__(cls
,
2205 def inner_product(self
, x
, y
):
2206 return _matrix_ip(x
,y
)
2209 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2211 The rank-n simple EJA consisting of complex Hermitian n-by-n
2212 matrices over the real numbers, the usual symmetric Jordan product,
2213 and the real-part-of-trace inner product. It has dimension `n^2` over
2218 The degree of this algebra is `n^2`::
2220 sage: set_random_seed()
2221 sage: n = ZZ.random_element(1,5)
2222 sage: J = ComplexHermitianEJA(n)
2223 sage: J.degree() == n^2
2226 The Jordan multiplication is what we think it is::
2228 sage: set_random_seed()
2229 sage: n = ZZ.random_element(1,5)
2230 sage: J = ComplexHermitianEJA(n)
2231 sage: x = J.random_element()
2232 sage: y = J.random_element()
2233 sage: actual = (x*y).natural_representation()
2234 sage: X = x.natural_representation()
2235 sage: Y = y.natural_representation()
2236 sage: expected = (X*Y + Y*X)/2
2237 sage: actual == expected
2239 sage: J(expected) == x*y
2244 def __classcall_private__(cls
, n
, field
=QQ
):
2245 S
= _complex_hermitian_basis(n
)
2246 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2248 fdeja
= super(ComplexHermitianEJA
, cls
)
2249 return fdeja
.__classcall
_private
__(cls
,
2255 def inner_product(self
, x
, y
):
2256 # Since a+bi on the diagonal is represented as
2261 # we'll double-count the "a" entries if we take the trace of
2263 return _matrix_ip(x
,y
)/2
2266 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2268 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2269 matrices, the usual symmetric Jordan product, and the
2270 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2275 The degree of this algebra is `n^2`::
2277 sage: set_random_seed()
2278 sage: n = ZZ.random_element(1,5)
2279 sage: J = QuaternionHermitianEJA(n)
2280 sage: J.degree() == 2*(n^2) - n
2283 The Jordan multiplication is what we think it is::
2285 sage: set_random_seed()
2286 sage: n = ZZ.random_element(1,5)
2287 sage: J = QuaternionHermitianEJA(n)
2288 sage: x = J.random_element()
2289 sage: y = J.random_element()
2290 sage: actual = (x*y).natural_representation()
2291 sage: X = x.natural_representation()
2292 sage: Y = y.natural_representation()
2293 sage: expected = (X*Y + Y*X)/2
2294 sage: actual == expected
2296 sage: J(expected) == x*y
2301 def __classcall_private__(cls
, n
, field
=QQ
):
2302 S
= _quaternion_hermitian_basis(n
)
2303 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2305 fdeja
= super(QuaternionHermitianEJA
, cls
)
2306 return fdeja
.__classcall
_private
__(cls
,
2312 def inner_product(self
, x
, y
):
2313 # Since a+bi+cj+dk on the diagonal is represented as
2315 # a + bi +cj + dk = [ a b c d]
2320 # we'll quadruple-count the "a" entries if we take the trace of
2322 return _matrix_ip(x
,y
)/4
2325 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2327 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2328 with the usual inner product and jordan product ``x*y =
2329 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2334 This multiplication table can be verified by hand::
2336 sage: J = JordanSpinEJA(4)
2337 sage: e0,e1,e2,e3 = J.gens()
2355 def __classcall_private__(cls
, n
, field
=QQ
):
2357 id_matrix
= identity_matrix(field
, n
)
2359 ei
= id_matrix
.column(i
)
2360 Qi
= zero_matrix(field
, n
)
2362 Qi
.set_column(0, ei
)
2363 Qi
+= diagonal_matrix(n
, [ei
[0]]*n
)
2364 # The addition of the diagonal matrix adds an extra ei[0] in the
2365 # upper-left corner of the matrix.
2366 Qi
[0,0] = Qi
[0,0] * ~
field(2)
2369 # The rank of the spin algebra is two, unless we're in a
2370 # one-dimensional ambient space (because the rank is bounded by
2371 # the ambient dimension).
2372 fdeja
= super(JordanSpinEJA
, cls
)
2373 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=min(n
,2))
2375 def inner_product(self
, x
, y
):
2376 return _usual_ip(x
,y
)