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.
10 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra
import FiniteDimensionalAlgebra
11 from sage
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra_element
import FiniteDimensionalAlgebraElement
12 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
13 from sage
.categories
.finite_dimensional_algebras_with_basis
import FiniteDimensionalAlgebrasWithBasis
14 from sage
.functions
.other
import sqrt
15 from sage
.matrix
.constructor
import matrix
16 from sage
.misc
.cachefunc
import cached_method
17 from sage
.misc
.prandom
import choice
18 from sage
.modules
.free_module
import VectorSpace
19 from sage
.modules
.free_module_element
import vector
20 from sage
.rings
.integer_ring
import ZZ
21 from sage
.rings
.number_field
.number_field
import QuadraticField
22 from sage
.rings
.polynomial
.polynomial_ring_constructor
import PolynomialRing
23 from sage
.rings
.rational_field
import QQ
24 from sage
.structure
.element
import is_Matrix
25 from sage
.structure
.category_object
import normalize_names
27 from mjo
.eja
.eja_operator
import FiniteDimensionalEuclideanJordanAlgebraOperator
30 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra
):
32 def __classcall_private__(cls
,
36 assume_associative
=False,
41 mult_table
= [b
.base_extend(field
) for b
in mult_table
]
44 if not (is_Matrix(b
) and b
.dimensions() == (n
, n
)):
45 raise ValueError("input is not a multiplication table")
46 mult_table
= tuple(mult_table
)
48 cat
= FiniteDimensionalAlgebrasWithBasis(field
)
49 cat
.or_subcategory(category
)
50 if assume_associative
:
51 cat
= cat
.Associative()
53 names
= normalize_names(n
, names
)
55 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, cls
)
56 return fda
.__classcall
__(cls
,
59 assume_associative
=assume_associative
,
63 natural_basis
=natural_basis
)
70 assume_associative
=False,
77 sage: from mjo.eja.eja_algebra import random_eja
81 By definition, Jordan multiplication commutes::
83 sage: set_random_seed()
84 sage: J = random_eja()
85 sage: x = J.random_element()
86 sage: y = J.random_element()
92 self
._natural
_basis
= natural_basis
93 self
._multiplication
_table
= mult_table
94 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
103 Return a string representation of ``self``.
105 fmt
= "Euclidean Jordan algebra of degree {} over {}"
106 return fmt
.format(self
.degree(), self
.base_ring())
109 def _a_regular_element(self
):
111 Guess a regular element. Needed to compute the basis for our
112 characteristic polynomial coefficients.
115 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
116 if not z
.is_regular():
117 raise ValueError("don't know a regular element")
122 def _charpoly_basis_space(self
):
124 Return the vector space spanned by the basis used in our
125 characteristic polynomial coefficients. This is used not only to
126 compute those coefficients, but also any time we need to
127 evaluate the coefficients (like when we compute the trace or
130 z
= self
._a
_regular
_element
()
131 V
= self
.vector_space()
132 V1
= V
.span_of_basis( (z
**k
).vector() for k
in range(self
.rank()) )
133 b
= (V1
.basis() + V1
.complement().basis())
134 return V
.span_of_basis(b
)
138 def _charpoly_coeff(self
, i
):
140 Return the coefficient polynomial "a_{i}" of this algebra's
141 general characteristic polynomial.
143 Having this be a separate cached method lets us compute and
144 store the trace/determinant (a_{r-1} and a_{0} respectively)
145 separate from the entire characteristic polynomial.
147 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
148 R
= A_of_x
.base_ring()
150 # Guaranteed by theory
153 # Danger: the in-place modification is done for performance
154 # reasons (reconstructing a matrix with huge polynomial
155 # entries is slow), but I don't know how cached_method works,
156 # so it's highly possible that we're modifying some global
157 # list variable by reference, here. In other words, you
158 # probably shouldn't call this method twice on the same
159 # algebra, at the same time, in two threads
160 Ai_orig
= A_of_x
.column(i
)
161 A_of_x
.set_column(i
,xr
)
162 numerator
= A_of_x
.det()
163 A_of_x
.set_column(i
,Ai_orig
)
165 # We're relying on the theory here to ensure that each a_i is
166 # indeed back in R, and the added negative signs are to make
167 # the whole charpoly expression sum to zero.
168 return R(-numerator
/detA
)
172 def _charpoly_matrix_system(self
):
174 Compute the matrix whose entries A_ij are polynomials in
175 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
176 corresponding to `x^r` and the determinent of the matrix A =
177 [A_ij]. In other words, all of the fixed (cachable) data needed
178 to compute the coefficients of the characteristic polynomial.
183 # Construct a new algebra over a multivariate polynomial ring...
184 names
= ['X' + str(i
) for i
in range(1,n
+1)]
185 R
= PolynomialRing(self
.base_ring(), names
)
186 J
= FiniteDimensionalEuclideanJordanAlgebra(R
,
187 self
._multiplication
_table
,
190 idmat
= matrix
.identity(J
.base_ring(), n
)
192 W
= self
._charpoly
_basis
_space
()
193 W
= W
.change_ring(R
.fraction_field())
195 # Starting with the standard coordinates x = (X1,X2,...,Xn)
196 # and then converting the entries to W-coordinates allows us
197 # to pass in the standard coordinates to the charpoly and get
198 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
201 # W.coordinates(x^2) eval'd at (standard z-coords)
205 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
207 # We want the middle equivalent thing in our matrix, but use
208 # the first equivalent thing instead so that we can pass in
209 # standard coordinates.
211 l1
= [matrix
.column(W
.coordinates((x
**k
).vector())) for k
in range(r
)]
212 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
213 A_of_x
= matrix
.block(R
, 1, n
, (l1
+ l2
))
214 xr
= W
.coordinates((x
**r
).vector())
215 return (A_of_x
, x
, xr
, A_of_x
.det())
219 def characteristic_polynomial(self
):
224 This implementation doesn't guarantee that the polynomial
225 denominator in the coefficients is not identically zero, so
226 theoretically it could crash. The way that this is handled
227 in e.g. Faraut and Koranyi is to use a basis that guarantees
228 the denominator is non-zero. But, doing so requires knowledge
229 of at least one regular element, and we don't even know how
230 to do that. The trade-off is that, if we use the standard basis,
231 the resulting polynomial will accept the "usual" coordinates. In
232 other words, we don't have to do a change of basis before e.g.
233 computing the trace or determinant.
237 sage: from mjo.eja.eja_algebra import JordanSpinEJA
241 The characteristic polynomial in the spin algebra is given in
242 Alizadeh, Example 11.11::
244 sage: J = JordanSpinEJA(3)
245 sage: p = J.characteristic_polynomial(); p
246 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
247 sage: xvec = J.one().vector()
255 # The list of coefficient polynomials a_1, a_2, ..., a_n.
256 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
258 # We go to a bit of trouble here to reorder the
259 # indeterminates, so that it's easier to evaluate the
260 # characteristic polynomial at x's coordinates and get back
261 # something in terms of t, which is what we want.
263 S
= PolynomialRing(self
.base_ring(),'t')
265 S
= PolynomialRing(S
, R
.variable_names())
268 # Note: all entries past the rth should be zero. The
269 # coefficient of the highest power (x^r) is 1, but it doesn't
270 # appear in the solution vector which contains coefficients
271 # for the other powers (to make them sum to x^r).
273 a
[r
] = 1 # corresponds to x^r
275 # When the rank is equal to the dimension, trying to
276 # assign a[r] goes out-of-bounds.
277 a
.append(1) # corresponds to x^r
279 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
282 def inner_product(self
, x
, y
):
284 The inner product associated with this Euclidean Jordan algebra.
286 Defaults to the trace inner product, but can be overridden by
287 subclasses if they are sure that the necessary properties are
292 sage: from mjo.eja.eja_algebra import random_eja
296 The inner product must satisfy its axiom for this algebra to truly
297 be a Euclidean Jordan Algebra::
299 sage: set_random_seed()
300 sage: J = random_eja()
301 sage: x = J.random_element()
302 sage: y = J.random_element()
303 sage: z = J.random_element()
304 sage: (x*y).inner_product(z) == y.inner_product(x*z)
308 if (not x
in self
) or (not y
in self
):
309 raise TypeError("arguments must live in this algebra")
310 return x
.trace_inner_product(y
)
313 def natural_basis(self
):
315 Return a more-natural representation of this algebra's basis.
317 Every finite-dimensional Euclidean Jordan Algebra is a direct
318 sum of five simple algebras, four of which comprise Hermitian
319 matrices. This method returns the original "natural" basis
320 for our underlying vector space. (Typically, the natural basis
321 is used to construct the multiplication table in the first place.)
323 Note that this will always return a matrix. The standard basis
324 in `R^n` will be returned as `n`-by-`1` column matrices.
328 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
329 ....: RealSymmetricEJA)
333 sage: J = RealSymmetricEJA(2)
336 sage: J.natural_basis()
344 sage: J = JordanSpinEJA(2)
347 sage: J.natural_basis()
354 if self
._natural
_basis
is None:
355 return tuple( b
.vector().column() for b
in self
.basis() )
357 return self
._natural
_basis
362 Return the rank of this EJA.
364 if self
._rank
is None:
365 raise ValueError("no rank specified at genesis")
370 def vector_space(self
):
372 Return the vector space that underlies this algebra.
376 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
380 sage: J = RealSymmetricEJA(2)
381 sage: J.vector_space()
382 Vector space of dimension 3 over Rational Field
385 return self
.zero().vector().parent().ambient_vector_space()
388 class Element(FiniteDimensionalAlgebraElement
):
390 An element of a Euclidean Jordan algebra.
395 Oh man, I should not be doing this. This hides the "disabled"
396 methods ``left_matrix`` and ``matrix`` from introspection;
397 in particular it removes them from tab-completion.
399 return filter(lambda s
: s
not in ['left_matrix', 'matrix'],
400 dir(self
.__class
__) )
403 def __init__(self
, A
, elt
=None):
408 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
412 The identity in `S^n` is converted to the identity in the EJA::
414 sage: J = RealSymmetricEJA(3)
415 sage: I = matrix.identity(QQ,3)
416 sage: J(I) == J.one()
419 This skew-symmetric matrix can't be represented in the EJA::
421 sage: J = RealSymmetricEJA(3)
422 sage: A = matrix(QQ,3, lambda i,j: i-j)
424 Traceback (most recent call last):
426 ArithmeticError: vector is not in free module
429 # Goal: if we're given a matrix, and if it lives in our
430 # parent algebra's "natural ambient space," convert it
431 # into an algebra element.
433 # The catch is, we make a recursive call after converting
434 # the given matrix into a vector that lives in the algebra.
435 # This we need to try the parent class initializer first,
436 # to avoid recursing forever if we're given something that
437 # already fits into the algebra, but also happens to live
438 # in the parent's "natural ambient space" (this happens with
441 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, elt
)
443 natural_basis
= A
.natural_basis()
444 if elt
in natural_basis
[0].matrix_space():
445 # Thanks for nothing! Matrix spaces aren't vector
446 # spaces in Sage, so we have to figure out its
447 # natural-basis coordinates ourselves.
448 V
= VectorSpace(elt
.base_ring(), elt
.nrows()**2)
449 W
= V
.span( _mat2vec(s
) for s
in natural_basis
)
450 coords
= W
.coordinates(_mat2vec(elt
))
451 FiniteDimensionalAlgebraElement
.__init
__(self
, A
, coords
)
453 def __pow__(self
, n
):
455 Return ``self`` raised to the power ``n``.
457 Jordan algebras are always power-associative; see for
458 example Faraut and Koranyi, Proposition II.1.2 (ii).
462 We have to override this because our superclass uses row vectors
463 instead of column vectors! We, on the other hand, assume column
468 sage: from mjo.eja.eja_algebra import random_eja
472 sage: set_random_seed()
473 sage: x = random_eja().random_element()
474 sage: x.operator()(x) == (x^2)
477 A few examples of power-associativity::
479 sage: set_random_seed()
480 sage: x = random_eja().random_element()
481 sage: x*(x*x)*(x*x) == x^5
483 sage: (x*x)*(x*x*x) == x^5
486 We also know that powers operator-commute (Koecher, Chapter
489 sage: set_random_seed()
490 sage: x = random_eja().random_element()
491 sage: m = ZZ.random_element(0,10)
492 sage: n = ZZ.random_element(0,10)
493 sage: Lxm = (x^m).operator()
494 sage: Lxn = (x^n).operator()
495 sage: Lxm*Lxn == Lxn*Lxm
500 return self
.parent().one()
504 return (self
.operator()**(n
-1))(self
)
507 def apply_univariate_polynomial(self
, p
):
509 Apply the univariate polynomial ``p`` to this element.
511 A priori, SageMath won't allow us to apply a univariate
512 polynomial to an element of an EJA, because we don't know
513 that EJAs are rings (they are usually not associative). Of
514 course, we know that EJAs are power-associative, so the
515 operation is ultimately kosher. This function sidesteps
516 the CAS to get the answer we want and expect.
520 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
525 sage: R = PolynomialRing(QQ, 't')
527 sage: p = t^4 - t^3 + 5*t - 2
528 sage: J = RealCartesianProductEJA(5)
529 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
534 We should always get back an element of the algebra::
536 sage: set_random_seed()
537 sage: p = PolynomialRing(QQ, 't').random_element()
538 sage: J = random_eja()
539 sage: x = J.random_element()
540 sage: x.apply_univariate_polynomial(p) in J
544 if len(p
.variables()) > 1:
545 raise ValueError("not a univariate polynomial")
548 # Convert the coeficcients to the parent's base ring,
549 # because a priori they might live in an (unnecessarily)
550 # larger ring for which P.sum() would fail below.
551 cs
= [ R(c
) for c
in p
.coefficients(sparse
=False) ]
552 return P
.sum( cs
[k
]*(self
**k
) for k
in range(len(cs
)) )
555 def characteristic_polynomial(self
):
557 Return the characteristic polynomial of this element.
561 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
565 The rank of `R^3` is three, and the minimal polynomial of
566 the identity element is `(t-1)` from which it follows that
567 the characteristic polynomial should be `(t-1)^3`::
569 sage: J = RealCartesianProductEJA(3)
570 sage: J.one().characteristic_polynomial()
571 t^3 - 3*t^2 + 3*t - 1
573 Likewise, the characteristic of the zero element in the
574 rank-three algebra `R^{n}` should be `t^{3}`::
576 sage: J = RealCartesianProductEJA(3)
577 sage: J.zero().characteristic_polynomial()
580 The characteristic polynomial of an element should evaluate
581 to zero on that element::
583 sage: set_random_seed()
584 sage: x = RealCartesianProductEJA(3).random_element()
585 sage: p = x.characteristic_polynomial()
586 sage: x.apply_univariate_polynomial(p)
590 p
= self
.parent().characteristic_polynomial()
591 return p(*self
.vector())
594 def inner_product(self
, other
):
596 Return the parent algebra's inner product of myself and ``other``.
600 sage: from mjo.eja.eja_algebra import (
601 ....: ComplexHermitianEJA,
603 ....: QuaternionHermitianEJA,
604 ....: RealSymmetricEJA,
609 The inner product in the Jordan spin algebra is the usual
610 inner product on `R^n` (this example only works because the
611 basis for the Jordan algebra is the standard basis in `R^n`)::
613 sage: J = JordanSpinEJA(3)
614 sage: x = vector(QQ,[1,2,3])
615 sage: y = vector(QQ,[4,5,6])
616 sage: x.inner_product(y)
618 sage: J(x).inner_product(J(y))
621 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
622 multiplication is the usual matrix multiplication in `S^n`,
623 so the inner product of the identity matrix with itself
626 sage: J = RealSymmetricEJA(3)
627 sage: J.one().inner_product(J.one())
630 Likewise, the inner product on `C^n` is `<X,Y> =
631 Re(trace(X*Y))`, where we must necessarily take the real
632 part because the product of Hermitian matrices may not be
635 sage: J = ComplexHermitianEJA(3)
636 sage: J.one().inner_product(J.one())
639 Ditto for the quaternions::
641 sage: J = QuaternionHermitianEJA(3)
642 sage: J.one().inner_product(J.one())
647 Ensure that we can always compute an inner product, and that
648 it gives us back a real number::
650 sage: set_random_seed()
651 sage: J = random_eja()
652 sage: x = J.random_element()
653 sage: y = J.random_element()
654 sage: x.inner_product(y) in RR
660 raise TypeError("'other' must live in the same algebra")
662 return P
.inner_product(self
, other
)
665 def operator_commutes_with(self
, other
):
667 Return whether or not this element operator-commutes
672 sage: from mjo.eja.eja_algebra import random_eja
676 The definition of a Jordan algebra says that any element
677 operator-commutes with its square::
679 sage: set_random_seed()
680 sage: x = random_eja().random_element()
681 sage: x.operator_commutes_with(x^2)
686 Test Lemma 1 from Chapter III of Koecher::
688 sage: set_random_seed()
689 sage: J = random_eja()
690 sage: u = J.random_element()
691 sage: v = J.random_element()
692 sage: lhs = u.operator_commutes_with(u*v)
693 sage: rhs = v.operator_commutes_with(u^2)
697 Test the first polarization identity from my notes, Koecher Chapter
698 III, or from Baes (2.3)::
700 sage: set_random_seed()
701 sage: J = random_eja()
702 sage: x = J.random_element()
703 sage: y = J.random_element()
704 sage: Lx = x.operator()
705 sage: Ly = y.operator()
706 sage: Lxx = (x*x).operator()
707 sage: Lxy = (x*y).operator()
708 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
711 Test the second polarization identity from my notes or from
714 sage: set_random_seed()
715 sage: J = random_eja()
716 sage: x = J.random_element()
717 sage: y = J.random_element()
718 sage: z = J.random_element()
719 sage: Lx = x.operator()
720 sage: Ly = y.operator()
721 sage: Lz = z.operator()
722 sage: Lzy = (z*y).operator()
723 sage: Lxy = (x*y).operator()
724 sage: Lxz = (x*z).operator()
725 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
728 Test the third polarization identity from my notes or from
731 sage: set_random_seed()
732 sage: J = random_eja()
733 sage: u = J.random_element()
734 sage: y = J.random_element()
735 sage: z = J.random_element()
736 sage: Lu = u.operator()
737 sage: Ly = y.operator()
738 sage: Lz = z.operator()
739 sage: Lzy = (z*y).operator()
740 sage: Luy = (u*y).operator()
741 sage: Luz = (u*z).operator()
742 sage: Luyz = (u*(y*z)).operator()
743 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
744 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
745 sage: bool(lhs == rhs)
749 if not other
in self
.parent():
750 raise TypeError("'other' must live in the same algebra")
759 Return my determinant, the product of my eigenvalues.
763 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
768 sage: J = JordanSpinEJA(2)
769 sage: e0,e1 = J.gens()
770 sage: x = sum( J.gens() )
776 sage: J = JordanSpinEJA(3)
777 sage: e0,e1,e2 = J.gens()
778 sage: x = sum( J.gens() )
784 An element is invertible if and only if its determinant is
787 sage: set_random_seed()
788 sage: x = random_eja().random_element()
789 sage: x.is_invertible() == (x.det() != 0)
795 p
= P
._charpoly
_coeff
(0)
796 # The _charpoly_coeff function already adds the factor of
797 # -1 to ensure that _charpoly_coeff(0) is really what
798 # appears in front of t^{0} in the charpoly. However,
799 # we want (-1)^r times THAT for the determinant.
800 return ((-1)**r
)*p(*self
.vector())
805 Return the Jordan-multiplicative inverse of this element.
809 We appeal to the quadratic representation as in Koecher's
810 Theorem 12 in Chapter III, Section 5.
814 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
819 The inverse in the spin factor algebra is given in Alizadeh's
822 sage: set_random_seed()
823 sage: n = ZZ.random_element(1,10)
824 sage: J = JordanSpinEJA(n)
825 sage: x = J.random_element()
826 sage: while not x.is_invertible():
827 ....: x = J.random_element()
828 sage: x_vec = x.vector()
830 sage: x_bar = x_vec[1:]
831 sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
832 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
833 sage: x_inverse = coeff*inv_vec
834 sage: x.inverse() == J(x_inverse)
839 The identity element is its own inverse::
841 sage: set_random_seed()
842 sage: J = random_eja()
843 sage: J.one().inverse() == J.one()
846 If an element has an inverse, it acts like one::
848 sage: set_random_seed()
849 sage: J = random_eja()
850 sage: x = J.random_element()
851 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
854 The inverse of the inverse is what we started with::
856 sage: set_random_seed()
857 sage: J = random_eja()
858 sage: x = J.random_element()
859 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
862 The zero element is never invertible::
864 sage: set_random_seed()
865 sage: J = random_eja().zero().inverse()
866 Traceback (most recent call last):
868 ValueError: element is not invertible
871 if not self
.is_invertible():
872 raise ValueError("element is not invertible")
874 return (~self
.quadratic_representation())(self
)
877 def is_invertible(self
):
879 Return whether or not this element is invertible.
881 We can't use the superclass method because it relies on
882 the algebra being associative.
886 The usual way to do this is to check if the determinant is
887 zero, but we need the characteristic polynomial for the
888 determinant. The minimal polynomial is a lot easier to get,
889 so we use Corollary 2 in Chapter V of Koecher to check
890 whether or not the paren't algebra's zero element is a root
891 of this element's minimal polynomial.
895 sage: from mjo.eja.eja_algebra import random_eja
899 The identity element is always invertible::
901 sage: set_random_seed()
902 sage: J = random_eja()
903 sage: J.one().is_invertible()
906 The zero element is never invertible::
908 sage: set_random_seed()
909 sage: J = random_eja()
910 sage: J.zero().is_invertible()
914 zero
= self
.parent().zero()
915 p
= self
.minimal_polynomial()
916 return not (p(zero
) == zero
)
919 def is_nilpotent(self
):
921 Return whether or not some power of this element is zero.
923 The superclass method won't work unless we're in an
924 associative algebra, and we aren't. However, we generate
925 an assocoative subalgebra and we're nilpotent there if and
926 only if we're nilpotent here (probably).
930 sage: from mjo.eja.eja_algebra import random_eja
934 The identity element is never nilpotent::
936 sage: set_random_seed()
937 sage: random_eja().one().is_nilpotent()
940 The additive identity is always nilpotent::
942 sage: set_random_seed()
943 sage: random_eja().zero().is_nilpotent()
947 # The element we're going to call "is_nilpotent()" on.
948 # Either myself, interpreted as an element of a finite-
949 # dimensional algebra, or an element of an associative
953 if self
.parent().is_associative():
954 elt
= FiniteDimensionalAlgebraElement(self
.parent(), self
)
956 V
= self
.span_of_powers()
957 assoc_subalg
= self
.subalgebra_generated_by()
958 # Mis-design warning: the basis used for span_of_powers()
959 # and subalgebra_generated_by() must be the same, and in
961 elt
= assoc_subalg(V
.coordinates(self
.vector()))
963 # Recursive call, but should work since elt lives in an
964 # associative algebra.
965 return elt
.is_nilpotent()
968 def is_regular(self
):
970 Return whether or not this is a regular element.
974 sage: from mjo.eja.eja_algebra import JordanSpinEJA
978 The identity element always has degree one, but any element
979 linearly-independent from it is regular::
981 sage: J = JordanSpinEJA(5)
982 sage: J.one().is_regular()
984 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
985 sage: for x in J.gens():
986 ....: (J.one() + x).is_regular()
994 return self
.degree() == self
.parent().rank()
999 Compute the degree of this element the straightforward way
1000 according to the definition; by appending powers to a list
1001 and figuring out its dimension (that is, whether or not
1002 they're linearly dependent).
1006 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1010 sage: J = JordanSpinEJA(4)
1011 sage: J.one().degree()
1013 sage: e0,e1,e2,e3 = J.gens()
1014 sage: (e0 - e1).degree()
1017 In the spin factor algebra (of rank two), all elements that
1018 aren't multiples of the identity are regular::
1020 sage: set_random_seed()
1021 sage: n = ZZ.random_element(1,10)
1022 sage: J = JordanSpinEJA(n)
1023 sage: x = J.random_element()
1024 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
1028 return self
.span_of_powers().dimension()
1031 def left_matrix(self
):
1033 Our parent class defines ``left_matrix`` and ``matrix``
1034 methods whose names are misleading. We don't want them.
1036 raise NotImplementedError("use operator().matrix() instead")
1038 matrix
= left_matrix
1041 def minimal_polynomial(self
):
1043 Return the minimal polynomial of this element,
1044 as a function of the variable `t`.
1048 We restrict ourselves to the associative subalgebra
1049 generated by this element, and then return the minimal
1050 polynomial of this element's operator matrix (in that
1051 subalgebra). This works by Baes Proposition 2.3.16.
1055 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1060 The minimal polynomial of the identity and zero elements are
1063 sage: set_random_seed()
1064 sage: J = random_eja()
1065 sage: J.one().minimal_polynomial()
1067 sage: J.zero().minimal_polynomial()
1070 The degree of an element is (by one definition) the degree
1071 of its minimal polynomial::
1073 sage: set_random_seed()
1074 sage: x = random_eja().random_element()
1075 sage: x.degree() == x.minimal_polynomial().degree()
1078 The minimal polynomial and the characteristic polynomial coincide
1079 and are known (see Alizadeh, Example 11.11) for all elements of
1080 the spin factor algebra that aren't scalar multiples of the
1083 sage: set_random_seed()
1084 sage: n = ZZ.random_element(2,10)
1085 sage: J = JordanSpinEJA(n)
1086 sage: y = J.random_element()
1087 sage: while y == y.coefficient(0)*J.one():
1088 ....: y = J.random_element()
1089 sage: y0 = y.vector()[0]
1090 sage: y_bar = y.vector()[1:]
1091 sage: actual = y.minimal_polynomial()
1092 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1093 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1094 sage: bool(actual == expected)
1097 The minimal polynomial should always kill its element::
1099 sage: set_random_seed()
1100 sage: x = random_eja().random_element()
1101 sage: p = x.minimal_polynomial()
1102 sage: x.apply_univariate_polynomial(p)
1106 V
= self
.span_of_powers()
1107 assoc_subalg
= self
.subalgebra_generated_by()
1108 # Mis-design warning: the basis used for span_of_powers()
1109 # and subalgebra_generated_by() must be the same, and in
1111 elt
= assoc_subalg(V
.coordinates(self
.vector()))
1112 return elt
.operator().minimal_polynomial()
1116 def natural_representation(self
):
1118 Return a more-natural representation of this element.
1120 Every finite-dimensional Euclidean Jordan Algebra is a
1121 direct sum of five simple algebras, four of which comprise
1122 Hermitian matrices. This method returns the original
1123 "natural" representation of this element as a Hermitian
1124 matrix, if it has one. If not, you get the usual representation.
1128 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1129 ....: QuaternionHermitianEJA)
1133 sage: J = ComplexHermitianEJA(3)
1136 sage: J.one().natural_representation()
1146 sage: J = QuaternionHermitianEJA(3)
1149 sage: J.one().natural_representation()
1150 [1 0 0 0 0 0 0 0 0 0 0 0]
1151 [0 1 0 0 0 0 0 0 0 0 0 0]
1152 [0 0 1 0 0 0 0 0 0 0 0 0]
1153 [0 0 0 1 0 0 0 0 0 0 0 0]
1154 [0 0 0 0 1 0 0 0 0 0 0 0]
1155 [0 0 0 0 0 1 0 0 0 0 0 0]
1156 [0 0 0 0 0 0 1 0 0 0 0 0]
1157 [0 0 0 0 0 0 0 1 0 0 0 0]
1158 [0 0 0 0 0 0 0 0 1 0 0 0]
1159 [0 0 0 0 0 0 0 0 0 1 0 0]
1160 [0 0 0 0 0 0 0 0 0 0 1 0]
1161 [0 0 0 0 0 0 0 0 0 0 0 1]
1164 B
= self
.parent().natural_basis()
1165 W
= B
[0].matrix_space()
1166 return W
.linear_combination(zip(self
.vector(), B
))
1171 Return the left-multiplication-by-this-element
1172 operator on the ambient algebra.
1176 sage: from mjo.eja.eja_algebra import random_eja
1180 sage: set_random_seed()
1181 sage: J = random_eja()
1182 sage: x = J.random_element()
1183 sage: y = J.random_element()
1184 sage: x.operator()(y) == x*y
1186 sage: y.operator()(x) == x*y
1191 fda_elt
= FiniteDimensionalAlgebraElement(P
, self
)
1192 return FiniteDimensionalEuclideanJordanAlgebraOperator(
1195 fda_elt
.matrix().transpose() )
1198 def quadratic_representation(self
, other
=None):
1200 Return the quadratic representation of this element.
1204 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1209 The explicit form in the spin factor algebra is given by
1210 Alizadeh's Example 11.12::
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_vec = x.vector()
1218 sage: x_bar = x_vec[1:]
1219 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1220 sage: B = 2*x0*x_bar.row()
1221 sage: C = 2*x0*x_bar.column()
1222 sage: D = matrix.identity(QQ, n-1)
1223 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1224 sage: D = D + 2*x_bar.tensor_product(x_bar)
1225 sage: Q = matrix.block(2,2,[A,B,C,D])
1226 sage: Q == x.quadratic_representation().matrix()
1229 Test all of the properties from Theorem 11.2 in Alizadeh::
1231 sage: set_random_seed()
1232 sage: J = random_eja()
1233 sage: x = J.random_element()
1234 sage: y = J.random_element()
1235 sage: Lx = x.operator()
1236 sage: Lxx = (x*x).operator()
1237 sage: Qx = x.quadratic_representation()
1238 sage: Qy = y.quadratic_representation()
1239 sage: Qxy = x.quadratic_representation(y)
1240 sage: Qex = J.one().quadratic_representation(x)
1241 sage: n = ZZ.random_element(10)
1242 sage: Qxn = (x^n).quadratic_representation()
1246 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1249 Property 2 (multiply on the right for :trac:`28272`):
1251 sage: alpha = QQ.random_element()
1252 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1257 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1260 sage: not x.is_invertible() or (
1263 ....: x.inverse().quadratic_representation() )
1266 sage: Qxy(J.one()) == x*y
1271 sage: not x.is_invertible() or (
1272 ....: x.quadratic_representation(x.inverse())*Qx
1273 ....: == Qx*x.quadratic_representation(x.inverse()) )
1276 sage: not x.is_invertible() or (
1277 ....: x.quadratic_representation(x.inverse())*Qx
1279 ....: 2*x.operator()*Qex - Qx )
1282 sage: 2*x.operator()*Qex - Qx == Lxx
1287 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1297 sage: not x.is_invertible() or (
1298 ....: Qx*x.inverse().operator() == Lx )
1303 sage: not x.operator_commutes_with(y) or (
1304 ....: Qx(y)^n == Qxn(y^n) )
1310 elif not other
in self
.parent():
1311 raise TypeError("'other' must live in the same algebra")
1314 M
= other
.operator()
1315 return ( L
*M
+ M
*L
- (self
*other
).operator() )
1318 def span_of_powers(self
):
1320 Return the vector space spanned by successive powers of
1323 # The dimension of the subalgebra can't be greater than
1324 # the big algebra, so just put everything into a list
1325 # and let span() get rid of the excess.
1327 # We do the extra ambient_vector_space() in case we're messing
1328 # with polynomials and the direct parent is a module.
1329 V
= self
.parent().vector_space()
1330 return V
.span( (self
**d
).vector() for d
in xrange(V
.dimension()) )
1333 def subalgebra_generated_by(self
):
1335 Return the associative subalgebra of the parent EJA generated
1340 sage: from mjo.eja.eja_algebra import random_eja
1344 sage: set_random_seed()
1345 sage: x = random_eja().random_element()
1346 sage: x.subalgebra_generated_by().is_associative()
1349 Squaring in the subalgebra should work the same as in
1352 sage: set_random_seed()
1353 sage: x = random_eja().random_element()
1354 sage: u = x.subalgebra_generated_by().random_element()
1355 sage: u.operator()(u) == u^2
1359 # First get the subspace spanned by the powers of myself...
1360 V
= self
.span_of_powers()
1361 F
= self
.base_ring()
1363 # Now figure out the entries of the right-multiplication
1364 # matrix for the successive basis elements b0, b1,... of
1367 for b_right
in V
.basis():
1368 eja_b_right
= self
.parent()(b_right
)
1370 # The first row of the right-multiplication matrix by
1371 # b1 is what we get if we apply that matrix to b1. The
1372 # second row of the right multiplication matrix by b1
1373 # is what we get when we apply that matrix to b2...
1375 # IMPORTANT: this assumes that all vectors are COLUMN
1376 # vectors, unlike our superclass (which uses row vectors).
1377 for b_left
in V
.basis():
1378 eja_b_left
= self
.parent()(b_left
)
1379 # Multiply in the original EJA, but then get the
1380 # coordinates from the subalgebra in terms of its
1382 this_row
= V
.coordinates((eja_b_left
*eja_b_right
).vector())
1383 b_right_rows
.append(this_row
)
1384 b_right_matrix
= matrix(F
, b_right_rows
)
1385 mats
.append(b_right_matrix
)
1387 # It's an algebra of polynomials in one element, and EJAs
1388 # are power-associative.
1390 # TODO: choose generator names intelligently.
1391 return FiniteDimensionalEuclideanJordanAlgebra(F
, mats
, assume_associative
=True, names
='f')
1394 def subalgebra_idempotent(self
):
1396 Find an idempotent in the associative subalgebra I generate
1397 using Proposition 2.3.5 in Baes.
1401 sage: from mjo.eja.eja_algebra import random_eja
1405 sage: set_random_seed()
1406 sage: J = random_eja()
1407 sage: x = J.random_element()
1408 sage: while x.is_nilpotent():
1409 ....: x = J.random_element()
1410 sage: c = x.subalgebra_idempotent()
1415 if self
.is_nilpotent():
1416 raise ValueError("this only works with non-nilpotent elements!")
1418 V
= self
.span_of_powers()
1419 J
= self
.subalgebra_generated_by()
1420 # Mis-design warning: the basis used for span_of_powers()
1421 # and subalgebra_generated_by() must be the same, and in
1423 u
= J(V
.coordinates(self
.vector()))
1425 # The image of the matrix of left-u^m-multiplication
1426 # will be minimal for some natural number s...
1428 minimal_dim
= V
.dimension()
1429 for i
in xrange(1, V
.dimension()):
1430 this_dim
= (u
**i
).operator().matrix().image().dimension()
1431 if this_dim
< minimal_dim
:
1432 minimal_dim
= this_dim
1435 # Now minimal_matrix should correspond to the smallest
1436 # non-zero subspace in Baes's (or really, Koecher's)
1439 # However, we need to restrict the matrix to work on the
1440 # subspace... or do we? Can't we just solve, knowing that
1441 # A(c) = u^(s+1) should have a solution in the big space,
1444 # Beware, solve_right() means that we're using COLUMN vectors.
1445 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1447 A
= u_next
.operator().matrix()
1448 c_coordinates
= A
.solve_right(u_next
.vector())
1450 # Now c_coordinates is the idempotent we want, but it's in
1451 # the coordinate system of the subalgebra.
1453 # We need the basis for J, but as elements of the parent algebra.
1455 basis
= [self
.parent(v
) for v
in V
.basis()]
1456 return self
.parent().linear_combination(zip(c_coordinates
, basis
))
1461 Return my trace, the sum of my eigenvalues.
1465 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1466 ....: RealCartesianProductEJA,
1471 sage: J = JordanSpinEJA(3)
1472 sage: x = sum(J.gens())
1478 sage: J = RealCartesianProductEJA(5)
1479 sage: J.one().trace()
1484 The trace of an element is a real number::
1486 sage: set_random_seed()
1487 sage: J = random_eja()
1488 sage: J.random_element().trace() in J.base_ring()
1494 p
= P
._charpoly
_coeff
(r
-1)
1495 # The _charpoly_coeff function already adds the factor of
1496 # -1 to ensure that _charpoly_coeff(r-1) is really what
1497 # appears in front of t^{r-1} in the charpoly. However,
1498 # we want the negative of THAT for the trace.
1499 return -p(*self
.vector())
1502 def trace_inner_product(self
, other
):
1504 Return the trace inner product of myself and ``other``.
1508 sage: from mjo.eja.eja_algebra import random_eja
1512 The trace inner product is commutative::
1514 sage: set_random_seed()
1515 sage: J = random_eja()
1516 sage: x = J.random_element(); y = J.random_element()
1517 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1520 The trace inner product is bilinear::
1522 sage: set_random_seed()
1523 sage: J = random_eja()
1524 sage: x = J.random_element()
1525 sage: y = J.random_element()
1526 sage: z = J.random_element()
1527 sage: a = QQ.random_element();
1528 sage: actual = (a*(x+z)).trace_inner_product(y)
1529 sage: expected = ( a*x.trace_inner_product(y) +
1530 ....: a*z.trace_inner_product(y) )
1531 sage: actual == expected
1533 sage: actual = x.trace_inner_product(a*(y+z))
1534 sage: expected = ( a*x.trace_inner_product(y) +
1535 ....: a*x.trace_inner_product(z) )
1536 sage: actual == expected
1539 The trace inner product satisfies the compatibility
1540 condition in the definition of a Euclidean Jordan algebra::
1542 sage: set_random_seed()
1543 sage: J = random_eja()
1544 sage: x = J.random_element()
1545 sage: y = J.random_element()
1546 sage: z = J.random_element()
1547 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1551 if not other
in self
.parent():
1552 raise TypeError("'other' must live in the same algebra")
1554 return (self
*other
).trace()
1557 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1559 Return the Euclidean Jordan Algebra corresponding to the set
1560 `R^n` under the Hadamard product.
1562 Note: this is nothing more than the Cartesian product of ``n``
1563 copies of the spin algebra. Once Cartesian product algebras
1564 are implemented, this can go.
1568 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
1572 This multiplication table can be verified by hand::
1574 sage: J = RealCartesianProductEJA(3)
1575 sage: e0,e1,e2 = J.gens()
1591 def __classcall_private__(cls
, n
, field
=QQ
):
1592 # The FiniteDimensionalAlgebra constructor takes a list of
1593 # matrices, the ith representing right multiplication by the ith
1594 # basis element in the vector space. So if e_1 = (1,0,0), then
1595 # right (Hadamard) multiplication of x by e_1 picks out the first
1596 # component of x; and likewise for the ith basis element e_i.
1597 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
1598 for i
in xrange(n
) ]
1600 fdeja
= super(RealCartesianProductEJA
, cls
)
1601 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=n
)
1603 def inner_product(self
, x
, y
):
1604 return _usual_ip(x
,y
)
1609 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1613 For now, we choose a random natural number ``n`` (greater than zero)
1614 and then give you back one of the following:
1616 * The cartesian product of the rational numbers ``n`` times; this is
1617 ``QQ^n`` with the Hadamard product.
1619 * The Jordan spin algebra on ``QQ^n``.
1621 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1624 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1625 in the space of ``2n``-by-``2n`` real symmetric matrices.
1627 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1628 in the space of ``4n``-by-``4n`` real symmetric matrices.
1630 Later this might be extended to return Cartesian products of the
1635 sage: from mjo.eja.eja_algebra import random_eja
1640 Euclidean Jordan algebra of degree...
1644 # The max_n component lets us choose different upper bounds on the
1645 # value "n" that gets passed to the constructor. This is needed
1646 # because e.g. R^{10} is reasonable to test, while the Hermitian
1647 # 10-by-10 quaternion matrices are not.
1648 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
1650 (RealSymmetricEJA
, 5),
1651 (ComplexHermitianEJA
, 4),
1652 (QuaternionHermitianEJA
, 3)])
1653 n
= ZZ
.random_element(1, max_n
)
1654 return constructor(n
, field
=QQ
)
1658 def _real_symmetric_basis(n
, field
=QQ
):
1660 Return a basis for the space of real symmetric n-by-n matrices.
1662 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1666 for j
in xrange(i
+1):
1667 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1671 # Beware, orthogonal but not normalized!
1672 Sij
= Eij
+ Eij
.transpose()
1677 def _complex_hermitian_basis(n
, field
=QQ
):
1679 Returns a basis for the space of complex Hermitian n-by-n matrices.
1683 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
1687 sage: set_random_seed()
1688 sage: n = ZZ.random_element(1,5)
1689 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1693 F
= QuadraticField(-1, 'I')
1696 # This is like the symmetric case, but we need to be careful:
1698 # * We want conjugate-symmetry, not just symmetry.
1699 # * The diagonal will (as a result) be real.
1703 for j
in xrange(i
+1):
1704 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1706 Sij
= _embed_complex_matrix(Eij
)
1709 # Beware, orthogonal but not normalized! The second one
1710 # has a minus because it's conjugated.
1711 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
1713 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
1718 def _quaternion_hermitian_basis(n
, field
=QQ
):
1720 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1724 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
1728 sage: set_random_seed()
1729 sage: n = ZZ.random_element(1,5)
1730 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1734 Q
= QuaternionAlgebra(QQ
,-1,-1)
1737 # This is like the symmetric case, but we need to be careful:
1739 # * We want conjugate-symmetry, not just symmetry.
1740 # * The diagonal will (as a result) be real.
1744 for j
in xrange(i
+1):
1745 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1747 Sij
= _embed_quaternion_matrix(Eij
)
1750 # Beware, orthogonal but not normalized! The second,
1751 # third, and fourth ones have a minus because they're
1753 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
1755 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
1757 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
1759 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
1765 return vector(m
.base_ring(), m
.list())
1768 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
1770 def _multiplication_table_from_matrix_basis(basis
):
1772 At least three of the five simple Euclidean Jordan algebras have the
1773 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1774 multiplication on the right is matrix multiplication. Given a basis
1775 for the underlying matrix space, this function returns a
1776 multiplication table (obtained by looping through the basis
1777 elements) for an algebra of those matrices. A reordered copy
1778 of the basis is also returned to work around the fact that
1779 the ``span()`` in this function will change the order of the basis
1780 from what we think it is, to... something else.
1782 # In S^2, for example, we nominally have four coordinates even
1783 # though the space is of dimension three only. The vector space V
1784 # is supposed to hold the entire long vector, and the subspace W
1785 # of V will be spanned by the vectors that arise from symmetric
1786 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1787 field
= basis
[0].base_ring()
1788 dimension
= basis
[0].nrows()
1790 V
= VectorSpace(field
, dimension
**2)
1791 W
= V
.span( _mat2vec(s
) for s
in basis
)
1793 # Taking the span above reorders our basis (thanks, jerk!) so we
1794 # need to put our "matrix basis" in the same order as the
1795 # (reordered) vector basis.
1796 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
1800 # Brute force the multiplication-by-s matrix by looping
1801 # through all elements of the basis and doing the computation
1802 # to find out what the corresponding row should be. BEWARE:
1803 # these multiplication tables won't be symmetric! It therefore
1804 # becomes REALLY IMPORTANT that the underlying algebra
1805 # constructor uses ROW vectors and not COLUMN vectors. That's
1806 # why we're computing rows here and not columns.
1809 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
1810 Q_rows
.append(W
.coordinates(this_row
))
1811 Q
= matrix(field
, W
.dimension(), Q_rows
)
1817 def _embed_complex_matrix(M
):
1819 Embed the n-by-n complex matrix ``M`` into the space of real
1820 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1821 bi` to the block matrix ``[[a,b],[-b,a]]``.
1825 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
1829 sage: F = QuadraticField(-1,'i')
1830 sage: x1 = F(4 - 2*i)
1831 sage: x2 = F(1 + 2*i)
1834 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1835 sage: _embed_complex_matrix(M)
1844 Embedding is a homomorphism (isomorphism, in fact)::
1846 sage: set_random_seed()
1847 sage: n = ZZ.random_element(5)
1848 sage: F = QuadraticField(-1, 'i')
1849 sage: X = random_matrix(F, n)
1850 sage: Y = random_matrix(F, n)
1851 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1852 sage: expected = _embed_complex_matrix(X*Y)
1853 sage: actual == expected
1859 raise ValueError("the matrix 'M' must be square")
1860 field
= M
.base_ring()
1865 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1867 # We can drop the imaginaries here.
1868 return matrix
.block(field
.base_ring(), n
, blocks
)
1871 def _unembed_complex_matrix(M
):
1873 The inverse of _embed_complex_matrix().
1877 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
1878 ....: _unembed_complex_matrix)
1882 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1883 ....: [-2, 1, -4, 3],
1884 ....: [ 9, 10, 11, 12],
1885 ....: [-10, 9, -12, 11] ])
1886 sage: _unembed_complex_matrix(A)
1888 [ 10*i + 9 12*i + 11]
1892 Unembedding is the inverse of embedding::
1894 sage: set_random_seed()
1895 sage: F = QuadraticField(-1, 'i')
1896 sage: M = random_matrix(F, 3)
1897 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1903 raise ValueError("the matrix 'M' must be square")
1904 if not n
.mod(2).is_zero():
1905 raise ValueError("the matrix 'M' must be a complex embedding")
1907 F
= QuadraticField(-1, 'i')
1910 # Go top-left to bottom-right (reading order), converting every
1911 # 2-by-2 block we see to a single complex element.
1913 for k
in xrange(n
/2):
1914 for j
in xrange(n
/2):
1915 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1916 if submat
[0,0] != submat
[1,1]:
1917 raise ValueError('bad on-diagonal submatrix')
1918 if submat
[0,1] != -submat
[1,0]:
1919 raise ValueError('bad off-diagonal submatrix')
1920 z
= submat
[0,0] + submat
[0,1]*i
1923 return matrix(F
, n
/2, elements
)
1926 def _embed_quaternion_matrix(M
):
1928 Embed the n-by-n quaternion matrix ``M`` into the space of real
1929 matrices of size 4n-by-4n by first sending each quaternion entry
1930 `z = a + bi + cj + dk` to the block-complex matrix
1931 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1936 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
1940 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1941 sage: i,j,k = Q.gens()
1942 sage: x = 1 + 2*i + 3*j + 4*k
1943 sage: M = matrix(Q, 1, [[x]])
1944 sage: _embed_quaternion_matrix(M)
1950 Embedding is a homomorphism (isomorphism, in fact)::
1952 sage: set_random_seed()
1953 sage: n = ZZ.random_element(5)
1954 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1955 sage: X = random_matrix(Q, n)
1956 sage: Y = random_matrix(Q, n)
1957 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1958 sage: expected = _embed_quaternion_matrix(X*Y)
1959 sage: actual == expected
1963 quaternions
= M
.base_ring()
1966 raise ValueError("the matrix 'M' must be square")
1968 F
= QuadraticField(-1, 'i')
1973 t
= z
.coefficient_tuple()
1978 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1979 [-c
+ d
*i
, a
- b
*i
]])
1980 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1982 # We should have real entries by now, so use the realest field
1983 # we've got for the return value.
1984 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1987 def _unembed_quaternion_matrix(M
):
1989 The inverse of _embed_quaternion_matrix().
1993 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1994 ....: _unembed_quaternion_matrix)
1998 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1999 ....: [-2, 1, -4, 3],
2000 ....: [-3, 4, 1, -2],
2001 ....: [-4, -3, 2, 1]])
2002 sage: _unembed_quaternion_matrix(M)
2003 [1 + 2*i + 3*j + 4*k]
2007 Unembedding is the inverse of embedding::
2009 sage: set_random_seed()
2010 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2011 sage: M = random_matrix(Q, 3)
2012 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
2018 raise ValueError("the matrix 'M' must be square")
2019 if not n
.mod(4).is_zero():
2020 raise ValueError("the matrix 'M' must be a complex embedding")
2022 Q
= QuaternionAlgebra(QQ
,-1,-1)
2025 # Go top-left to bottom-right (reading order), converting every
2026 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2029 for l
in xrange(n
/4):
2030 for m
in xrange(n
/4):
2031 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
2032 if submat
[0,0] != submat
[1,1].conjugate():
2033 raise ValueError('bad on-diagonal submatrix')
2034 if submat
[0,1] != -submat
[1,0].conjugate():
2035 raise ValueError('bad off-diagonal submatrix')
2036 z
= submat
[0,0].real() + submat
[0,0].imag()*i
2037 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
2040 return matrix(Q
, n
/4, elements
)
2043 # The usual inner product on R^n.
2045 return x
.vector().inner_product(y
.vector())
2047 # The inner product used for the real symmetric simple EJA.
2048 # We keep it as a separate function because e.g. the complex
2049 # algebra uses the same inner product, except divided by 2.
2050 def _matrix_ip(X
,Y
):
2051 X_mat
= X
.natural_representation()
2052 Y_mat
= Y
.natural_representation()
2053 return (X_mat
*Y_mat
).trace()
2056 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2058 The rank-n simple EJA consisting of real symmetric n-by-n
2059 matrices, the usual symmetric Jordan product, and the trace inner
2060 product. It has dimension `(n^2 + n)/2` over the reals.
2064 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
2068 sage: J = RealSymmetricEJA(2)
2069 sage: e0, e1, e2 = J.gens()
2079 The degree of this algebra is `(n^2 + n) / 2`::
2081 sage: set_random_seed()
2082 sage: n = ZZ.random_element(1,5)
2083 sage: J = RealSymmetricEJA(n)
2084 sage: J.degree() == (n^2 + n)/2
2087 The Jordan multiplication is what we think it is::
2089 sage: set_random_seed()
2090 sage: n = ZZ.random_element(1,5)
2091 sage: J = RealSymmetricEJA(n)
2092 sage: x = J.random_element()
2093 sage: y = J.random_element()
2094 sage: actual = (x*y).natural_representation()
2095 sage: X = x.natural_representation()
2096 sage: Y = y.natural_representation()
2097 sage: expected = (X*Y + Y*X)/2
2098 sage: actual == expected
2100 sage: J(expected) == x*y
2105 def __classcall_private__(cls
, n
, field
=QQ
):
2106 S
= _real_symmetric_basis(n
, field
=field
)
2107 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2109 fdeja
= super(RealSymmetricEJA
, cls
)
2110 return fdeja
.__classcall
_private
__(cls
,
2116 def inner_product(self
, x
, y
):
2117 return _matrix_ip(x
,y
)
2120 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2122 The rank-n simple EJA consisting of complex Hermitian n-by-n
2123 matrices over the real numbers, the usual symmetric Jordan product,
2124 and the real-part-of-trace inner product. It has dimension `n^2` over
2129 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2133 The degree of this algebra is `n^2`::
2135 sage: set_random_seed()
2136 sage: n = ZZ.random_element(1,5)
2137 sage: J = ComplexHermitianEJA(n)
2138 sage: J.degree() == n^2
2141 The Jordan multiplication is what we think it is::
2143 sage: set_random_seed()
2144 sage: n = ZZ.random_element(1,5)
2145 sage: J = ComplexHermitianEJA(n)
2146 sage: x = J.random_element()
2147 sage: y = J.random_element()
2148 sage: actual = (x*y).natural_representation()
2149 sage: X = x.natural_representation()
2150 sage: Y = y.natural_representation()
2151 sage: expected = (X*Y + Y*X)/2
2152 sage: actual == expected
2154 sage: J(expected) == x*y
2159 def __classcall_private__(cls
, n
, field
=QQ
):
2160 S
= _complex_hermitian_basis(n
)
2161 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2163 fdeja
= super(ComplexHermitianEJA
, cls
)
2164 return fdeja
.__classcall
_private
__(cls
,
2170 def inner_product(self
, x
, y
):
2171 # Since a+bi on the diagonal is represented as
2176 # we'll double-count the "a" entries if we take the trace of
2178 return _matrix_ip(x
,y
)/2
2181 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2183 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2184 matrices, the usual symmetric Jordan product, and the
2185 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2190 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2194 The degree of this algebra is `n^2`::
2196 sage: set_random_seed()
2197 sage: n = ZZ.random_element(1,5)
2198 sage: J = QuaternionHermitianEJA(n)
2199 sage: J.degree() == 2*(n^2) - n
2202 The Jordan multiplication is what we think it is::
2204 sage: set_random_seed()
2205 sage: n = ZZ.random_element(1,5)
2206 sage: J = QuaternionHermitianEJA(n)
2207 sage: x = J.random_element()
2208 sage: y = J.random_element()
2209 sage: actual = (x*y).natural_representation()
2210 sage: X = x.natural_representation()
2211 sage: Y = y.natural_representation()
2212 sage: expected = (X*Y + Y*X)/2
2213 sage: actual == expected
2215 sage: J(expected) == x*y
2220 def __classcall_private__(cls
, n
, field
=QQ
):
2221 S
= _quaternion_hermitian_basis(n
)
2222 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
2224 fdeja
= super(QuaternionHermitianEJA
, cls
)
2225 return fdeja
.__classcall
_private
__(cls
,
2231 def inner_product(self
, x
, y
):
2232 # Since a+bi+cj+dk on the diagonal is represented as
2234 # a + bi +cj + dk = [ a b c d]
2239 # we'll quadruple-count the "a" entries if we take the trace of
2241 return _matrix_ip(x
,y
)/4
2244 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2246 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2247 with the usual inner product and jordan product ``x*y =
2248 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2253 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2257 This multiplication table can be verified by hand::
2259 sage: J = JordanSpinEJA(4)
2260 sage: e0,e1,e2,e3 = J.gens()
2278 def __classcall_private__(cls
, n
, field
=QQ
):
2280 id_matrix
= matrix
.identity(field
, n
)
2282 ei
= id_matrix
.column(i
)
2283 Qi
= matrix
.zero(field
, n
)
2285 Qi
.set_column(0, ei
)
2286 Qi
+= matrix
.diagonal(n
, [ei
[0]]*n
)
2287 # The addition of the diagonal matrix adds an extra ei[0] in the
2288 # upper-left corner of the matrix.
2289 Qi
[0,0] = Qi
[0,0] * ~
field(2)
2292 # The rank of the spin algebra is two, unless we're in a
2293 # one-dimensional ambient space (because the rank is bounded by
2294 # the ambient dimension).
2295 fdeja
= super(JordanSpinEJA
, cls
)
2296 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=min(n
,2))
2298 def inner_product(self
, x
, y
):
2299 return _usual_ip(x
,y
)