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
.algebras
.finite_dimensional_algebras
.finite_dimensional_algebra
import FiniteDimensionalAlgebra
9 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
10 from sage
.categories
.finite_dimensional_algebras_with_basis
import FiniteDimensionalAlgebrasWithBasis
11 from sage
.matrix
.constructor
import matrix
12 from sage
.misc
.cachefunc
import cached_method
13 from sage
.misc
.prandom
import choice
14 from sage
.modules
.free_module
import VectorSpace
15 from sage
.rings
.integer_ring
import ZZ
16 from sage
.rings
.number_field
.number_field
import QuadraticField
17 from sage
.rings
.polynomial
.polynomial_ring_constructor
import PolynomialRing
18 from sage
.rings
.rational_field
import QQ
19 from sage
.structure
.element
import is_Matrix
20 from sage
.structure
.category_object
import normalize_names
22 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
23 from mjo
.eja
.eja_utils
import _vec2mat
, _mat2vec
25 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra
):
27 def __classcall_private__(cls
,
32 assume_associative
=False,
36 mult_table
= [b
.base_extend(field
) for b
in mult_table
]
39 if not (is_Matrix(b
) and b
.dimensions() == (n
, n
)):
40 raise ValueError("input is not a multiplication table")
41 mult_table
= tuple(mult_table
)
43 cat
= FiniteDimensionalAlgebrasWithBasis(field
)
44 cat
.or_subcategory(category
)
45 if assume_associative
:
46 cat
= cat
.Associative()
48 names
= normalize_names(n
, names
)
50 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, cls
)
51 return fda
.__classcall
__(cls
,
55 assume_associative
=assume_associative
,
58 natural_basis
=natural_basis
)
66 assume_associative
=False,
72 sage: from mjo.eja.eja_algebra import random_eja
76 By definition, Jordan multiplication commutes::
78 sage: set_random_seed()
79 sage: J = random_eja()
80 sage: x = J.random_element()
81 sage: y = J.random_element()
87 self
._natural
_basis
= natural_basis
88 self
._multiplication
_table
= mult_table
89 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
98 Return a string representation of ``self``.
102 sage: from mjo.eja.eja_algebra import JordanSpinEJA
106 Ensure that it says what we think it says::
108 sage: JordanSpinEJA(2, field=QQ)
109 Euclidean Jordan algebra of degree 2 over Rational Field
110 sage: JordanSpinEJA(3, field=RDF)
111 Euclidean Jordan algebra of degree 3 over Real Double Field
114 fmt
= "Euclidean Jordan algebra of degree {} over {}"
115 return fmt
.format(self
.degree(), self
.base_ring())
118 def _a_regular_element(self
):
120 Guess a regular element. Needed to compute the basis for our
121 characteristic polynomial coefficients.
125 sage: from mjo.eja.eja_algebra import random_eja
129 Ensure that this hacky method succeeds for every algebra that we
130 know how to construct::
132 sage: set_random_seed()
133 sage: J = random_eja()
134 sage: J._a_regular_element().is_regular()
139 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
140 if not z
.is_regular():
141 raise ValueError("don't know a regular element")
146 def _charpoly_basis_space(self
):
148 Return the vector space spanned by the basis used in our
149 characteristic polynomial coefficients. This is used not only to
150 compute those coefficients, but also any time we need to
151 evaluate the coefficients (like when we compute the trace or
154 z
= self
._a
_regular
_element
()
155 V
= self
.vector_space()
156 V1
= V
.span_of_basis( (z
**k
).vector() for k
in range(self
.rank()) )
157 b
= (V1
.basis() + V1
.complement().basis())
158 return V
.span_of_basis(b
)
162 def _charpoly_coeff(self
, i
):
164 Return the coefficient polynomial "a_{i}" of this algebra's
165 general characteristic polynomial.
167 Having this be a separate cached method lets us compute and
168 store the trace/determinant (a_{r-1} and a_{0} respectively)
169 separate from the entire characteristic polynomial.
171 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
172 R
= A_of_x
.base_ring()
174 # Guaranteed by theory
177 # Danger: the in-place modification is done for performance
178 # reasons (reconstructing a matrix with huge polynomial
179 # entries is slow), but I don't know how cached_method works,
180 # so it's highly possible that we're modifying some global
181 # list variable by reference, here. In other words, you
182 # probably shouldn't call this method twice on the same
183 # algebra, at the same time, in two threads
184 Ai_orig
= A_of_x
.column(i
)
185 A_of_x
.set_column(i
,xr
)
186 numerator
= A_of_x
.det()
187 A_of_x
.set_column(i
,Ai_orig
)
189 # We're relying on the theory here to ensure that each a_i is
190 # indeed back in R, and the added negative signs are to make
191 # the whole charpoly expression sum to zero.
192 return R(-numerator
/detA
)
196 def _charpoly_matrix_system(self
):
198 Compute the matrix whose entries A_ij are polynomials in
199 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
200 corresponding to `x^r` and the determinent of the matrix A =
201 [A_ij]. In other words, all of the fixed (cachable) data needed
202 to compute the coefficients of the characteristic polynomial.
207 # Construct a new algebra over a multivariate polynomial ring...
208 names
= ['X' + str(i
) for i
in range(1,n
+1)]
209 R
= PolynomialRing(self
.base_ring(), names
)
210 J
= FiniteDimensionalEuclideanJordanAlgebra(R
,
211 self
._multiplication
_table
,
214 idmat
= matrix
.identity(J
.base_ring(), n
)
216 W
= self
._charpoly
_basis
_space
()
217 W
= W
.change_ring(R
.fraction_field())
219 # Starting with the standard coordinates x = (X1,X2,...,Xn)
220 # and then converting the entries to W-coordinates allows us
221 # to pass in the standard coordinates to the charpoly and get
222 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
225 # W.coordinates(x^2) eval'd at (standard z-coords)
229 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
231 # We want the middle equivalent thing in our matrix, but use
232 # the first equivalent thing instead so that we can pass in
233 # standard coordinates.
236 # Handle the zeroth power separately, because computing
237 # the unit element in J is mathematically suspect.
238 x0
= W
.coordinates(self
.one().vector())
239 l1
= [ matrix
.column(x0
) ]
240 l1
+= [ matrix
.column(W
.coordinates((x
**k
).vector()))
241 for k
in range(1,r
) ]
242 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
243 A_of_x
= matrix
.block(R
, 1, n
, (l1
+ l2
))
244 xr
= W
.coordinates((x
**r
).vector())
245 return (A_of_x
, x
, xr
, A_of_x
.det())
249 def characteristic_polynomial(self
):
251 Return a characteristic polynomial that works for all elements
254 The resulting polynomial has `n+1` variables, where `n` is the
255 dimension of this algebra. The first `n` variables correspond to
256 the coordinates of an algebra element: when evaluated at the
257 coordinates of an algebra element with respect to a certain
258 basis, the result is a univariate polynomial (in the one
259 remaining variable ``t``), namely the characteristic polynomial
264 sage: from mjo.eja.eja_algebra import JordanSpinEJA
268 The characteristic polynomial in the spin algebra is given in
269 Alizadeh, Example 11.11::
271 sage: J = JordanSpinEJA(3)
272 sage: p = J.characteristic_polynomial(); p
273 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
274 sage: xvec = J.one().vector()
282 # The list of coefficient polynomials a_1, a_2, ..., a_n.
283 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
285 # We go to a bit of trouble here to reorder the
286 # indeterminates, so that it's easier to evaluate the
287 # characteristic polynomial at x's coordinates and get back
288 # something in terms of t, which is what we want.
290 S
= PolynomialRing(self
.base_ring(),'t')
292 S
= PolynomialRing(S
, R
.variable_names())
295 # Note: all entries past the rth should be zero. The
296 # coefficient of the highest power (x^r) is 1, but it doesn't
297 # appear in the solution vector which contains coefficients
298 # for the other powers (to make them sum to x^r).
300 a
[r
] = 1 # corresponds to x^r
302 # When the rank is equal to the dimension, trying to
303 # assign a[r] goes out-of-bounds.
304 a
.append(1) # corresponds to x^r
306 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
309 def inner_product(self
, x
, y
):
311 The inner product associated with this Euclidean Jordan algebra.
313 Defaults to the trace inner product, but can be overridden by
314 subclasses if they are sure that the necessary properties are
319 sage: from mjo.eja.eja_algebra import random_eja
323 The inner product must satisfy its axiom for this algebra to truly
324 be a Euclidean Jordan Algebra::
326 sage: set_random_seed()
327 sage: J = random_eja()
328 sage: x = J.random_element()
329 sage: y = J.random_element()
330 sage: z = J.random_element()
331 sage: (x*y).inner_product(z) == y.inner_product(x*z)
335 if (not x
in self
) or (not y
in self
):
336 raise TypeError("arguments must live in this algebra")
337 return x
.trace_inner_product(y
)
340 def natural_basis(self
):
342 Return a more-natural representation of this algebra's basis.
344 Every finite-dimensional Euclidean Jordan Algebra is a direct
345 sum of five simple algebras, four of which comprise Hermitian
346 matrices. This method returns the original "natural" basis
347 for our underlying vector space. (Typically, the natural basis
348 is used to construct the multiplication table in the first place.)
350 Note that this will always return a matrix. The standard basis
351 in `R^n` will be returned as `n`-by-`1` column matrices.
355 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
356 ....: RealSymmetricEJA)
360 sage: J = RealSymmetricEJA(2)
363 sage: J.natural_basis()
371 sage: J = JordanSpinEJA(2)
374 sage: J.natural_basis()
381 if self
._natural
_basis
is None:
382 return tuple( b
.vector().column() for b
in self
.basis() )
384 return self
._natural
_basis
390 Return the unit element of this algebra.
394 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
399 sage: J = RealCartesianProductEJA(5)
401 e0 + e1 + e2 + e3 + e4
405 The identity element acts like the identity::
407 sage: set_random_seed()
408 sage: J = random_eja()
409 sage: x = J.random_element()
410 sage: J.one()*x == x and x*J.one() == x
413 The matrix of the unit element's operator is the identity::
415 sage: set_random_seed()
416 sage: J = random_eja()
417 sage: actual = J.one().operator().matrix()
418 sage: expected = matrix.identity(J.base_ring(), J.dimension())
419 sage: actual == expected
423 # We can brute-force compute the matrices of the operators
424 # that correspond to the basis elements of this algebra.
425 # If some linear combination of those basis elements is the
426 # algebra identity, then the same linear combination of
427 # their matrices has to be the identity matrix.
429 # Of course, matrices aren't vectors in sage, so we have to
430 # appeal to the "long vectors" isometry.
431 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
433 # Now we use basis linear algebra to find the coefficients,
434 # of the matrices-as-vectors-linear-combination, which should
435 # work for the original algebra basis too.
436 A
= matrix
.column(self
.base_ring(), oper_vecs
)
438 # We used the isometry on the left-hand side already, but we
439 # still need to do it for the right-hand side. Recall that we
440 # wanted something that summed to the identity matrix.
441 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
443 # Now if there's an identity element in the algebra, this should work.
444 coeffs
= A
.solve_right(b
)
445 return self
.linear_combination(zip(coeffs
,self
.gens()))
450 Return the rank of this EJA.
454 The author knows of no algorithm to compute the rank of an EJA
455 where only the multiplication table is known. In lieu of one, we
456 require the rank to be specified when the algebra is created,
457 and simply pass along that number here.
461 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
462 ....: RealSymmetricEJA,
463 ....: ComplexHermitianEJA,
464 ....: QuaternionHermitianEJA,
469 The rank of the Jordan spin algebra is always two::
471 sage: JordanSpinEJA(2).rank()
473 sage: JordanSpinEJA(3).rank()
475 sage: JordanSpinEJA(4).rank()
478 The rank of the `n`-by-`n` Hermitian real, complex, or
479 quaternion matrices is `n`::
481 sage: RealSymmetricEJA(2).rank()
483 sage: ComplexHermitianEJA(2).rank()
485 sage: QuaternionHermitianEJA(2).rank()
487 sage: RealSymmetricEJA(5).rank()
489 sage: ComplexHermitianEJA(5).rank()
491 sage: QuaternionHermitianEJA(5).rank()
496 Ensure that every EJA that we know how to construct has a
497 positive integer rank::
499 sage: set_random_seed()
500 sage: r = random_eja().rank()
501 sage: r in ZZ and r > 0
508 def vector_space(self
):
510 Return the vector space that underlies this algebra.
514 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
518 sage: J = RealSymmetricEJA(2)
519 sage: J.vector_space()
520 Vector space of dimension 3 over Rational Field
523 return self
.zero().vector().parent().ambient_vector_space()
526 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
529 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
531 Return the Euclidean Jordan Algebra corresponding to the set
532 `R^n` under the Hadamard product.
534 Note: this is nothing more than the Cartesian product of ``n``
535 copies of the spin algebra. Once Cartesian product algebras
536 are implemented, this can go.
540 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
544 This multiplication table can be verified by hand::
546 sage: J = RealCartesianProductEJA(3)
547 sage: e0,e1,e2 = J.gens()
563 def __classcall_private__(cls
, n
, field
=QQ
):
564 # The FiniteDimensionalAlgebra constructor takes a list of
565 # matrices, the ith representing right multiplication by the ith
566 # basis element in the vector space. So if e_1 = (1,0,0), then
567 # right (Hadamard) multiplication of x by e_1 picks out the first
568 # component of x; and likewise for the ith basis element e_i.
569 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
572 fdeja
= super(RealCartesianProductEJA
, cls
)
573 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=n
)
575 def inner_product(self
, x
, y
):
576 return _usual_ip(x
,y
)
581 Return a "random" finite-dimensional Euclidean Jordan Algebra.
585 For now, we choose a random natural number ``n`` (greater than zero)
586 and then give you back one of the following:
588 * The cartesian product of the rational numbers ``n`` times; this is
589 ``QQ^n`` with the Hadamard product.
591 * The Jordan spin algebra on ``QQ^n``.
593 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
596 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
597 in the space of ``2n``-by-``2n`` real symmetric matrices.
599 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
600 in the space of ``4n``-by-``4n`` real symmetric matrices.
602 Later this might be extended to return Cartesian products of the
607 sage: from mjo.eja.eja_algebra import random_eja
612 Euclidean Jordan algebra of degree...
616 # The max_n component lets us choose different upper bounds on the
617 # value "n" that gets passed to the constructor. This is needed
618 # because e.g. R^{10} is reasonable to test, while the Hermitian
619 # 10-by-10 quaternion matrices are not.
620 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
622 (RealSymmetricEJA
, 5),
623 (ComplexHermitianEJA
, 4),
624 (QuaternionHermitianEJA
, 3)])
625 n
= ZZ
.random_element(1, max_n
)
626 return constructor(n
, field
=QQ
)
630 def _real_symmetric_basis(n
, field
=QQ
):
632 Return a basis for the space of real symmetric n-by-n matrices.
634 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
638 for j
in xrange(i
+1):
639 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
643 # Beware, orthogonal but not normalized!
644 Sij
= Eij
+ Eij
.transpose()
649 def _complex_hermitian_basis(n
, field
=QQ
):
651 Returns a basis for the space of complex Hermitian n-by-n matrices.
655 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
659 sage: set_random_seed()
660 sage: n = ZZ.random_element(1,5)
661 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
665 F
= QuadraticField(-1, 'I')
668 # This is like the symmetric case, but we need to be careful:
670 # * We want conjugate-symmetry, not just symmetry.
671 # * The diagonal will (as a result) be real.
675 for j
in xrange(i
+1):
676 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
678 Sij
= _embed_complex_matrix(Eij
)
681 # Beware, orthogonal but not normalized! The second one
682 # has a minus because it's conjugated.
683 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
685 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
690 def _quaternion_hermitian_basis(n
, field
=QQ
):
692 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
696 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
700 sage: set_random_seed()
701 sage: n = ZZ.random_element(1,5)
702 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
706 Q
= QuaternionAlgebra(QQ
,-1,-1)
709 # This is like the symmetric case, but we need to be careful:
711 # * We want conjugate-symmetry, not just symmetry.
712 # * The diagonal will (as a result) be real.
716 for j
in xrange(i
+1):
717 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
719 Sij
= _embed_quaternion_matrix(Eij
)
722 # Beware, orthogonal but not normalized! The second,
723 # third, and fourth ones have a minus because they're
725 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
727 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
729 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
731 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
737 def _multiplication_table_from_matrix_basis(basis
):
739 At least three of the five simple Euclidean Jordan algebras have the
740 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
741 multiplication on the right is matrix multiplication. Given a basis
742 for the underlying matrix space, this function returns a
743 multiplication table (obtained by looping through the basis
744 elements) for an algebra of those matrices. A reordered copy
745 of the basis is also returned to work around the fact that
746 the ``span()`` in this function will change the order of the basis
747 from what we think it is, to... something else.
749 # In S^2, for example, we nominally have four coordinates even
750 # though the space is of dimension three only. The vector space V
751 # is supposed to hold the entire long vector, and the subspace W
752 # of V will be spanned by the vectors that arise from symmetric
753 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
754 field
= basis
[0].base_ring()
755 dimension
= basis
[0].nrows()
757 V
= VectorSpace(field
, dimension
**2)
758 W
= V
.span( _mat2vec(s
) for s
in basis
)
760 # Taking the span above reorders our basis (thanks, jerk!) so we
761 # need to put our "matrix basis" in the same order as the
762 # (reordered) vector basis.
763 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
767 # Brute force the multiplication-by-s matrix by looping
768 # through all elements of the basis and doing the computation
769 # to find out what the corresponding row should be. BEWARE:
770 # these multiplication tables won't be symmetric! It therefore
771 # becomes REALLY IMPORTANT that the underlying algebra
772 # constructor uses ROW vectors and not COLUMN vectors. That's
773 # why we're computing rows here and not columns.
776 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
777 Q_rows
.append(W
.coordinates(this_row
))
778 Q
= matrix(field
, W
.dimension(), Q_rows
)
784 def _embed_complex_matrix(M
):
786 Embed the n-by-n complex matrix ``M`` into the space of real
787 matrices of size 2n-by-2n via the map the sends each entry `z = a +
788 bi` to the block matrix ``[[a,b],[-b,a]]``.
792 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
796 sage: F = QuadraticField(-1,'i')
797 sage: x1 = F(4 - 2*i)
798 sage: x2 = F(1 + 2*i)
801 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
802 sage: _embed_complex_matrix(M)
811 Embedding is a homomorphism (isomorphism, in fact)::
813 sage: set_random_seed()
814 sage: n = ZZ.random_element(5)
815 sage: F = QuadraticField(-1, 'i')
816 sage: X = random_matrix(F, n)
817 sage: Y = random_matrix(F, n)
818 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
819 sage: expected = _embed_complex_matrix(X*Y)
820 sage: actual == expected
826 raise ValueError("the matrix 'M' must be square")
827 field
= M
.base_ring()
832 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
834 # We can drop the imaginaries here.
835 return matrix
.block(field
.base_ring(), n
, blocks
)
838 def _unembed_complex_matrix(M
):
840 The inverse of _embed_complex_matrix().
844 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
845 ....: _unembed_complex_matrix)
849 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
850 ....: [-2, 1, -4, 3],
851 ....: [ 9, 10, 11, 12],
852 ....: [-10, 9, -12, 11] ])
853 sage: _unembed_complex_matrix(A)
855 [ 10*i + 9 12*i + 11]
859 Unembedding is the inverse of embedding::
861 sage: set_random_seed()
862 sage: F = QuadraticField(-1, 'i')
863 sage: M = random_matrix(F, 3)
864 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
870 raise ValueError("the matrix 'M' must be square")
871 if not n
.mod(2).is_zero():
872 raise ValueError("the matrix 'M' must be a complex embedding")
874 F
= QuadraticField(-1, 'i')
877 # Go top-left to bottom-right (reading order), converting every
878 # 2-by-2 block we see to a single complex element.
880 for k
in xrange(n
/2):
881 for j
in xrange(n
/2):
882 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
883 if submat
[0,0] != submat
[1,1]:
884 raise ValueError('bad on-diagonal submatrix')
885 if submat
[0,1] != -submat
[1,0]:
886 raise ValueError('bad off-diagonal submatrix')
887 z
= submat
[0,0] + submat
[0,1]*i
890 return matrix(F
, n
/2, elements
)
893 def _embed_quaternion_matrix(M
):
895 Embed the n-by-n quaternion matrix ``M`` into the space of real
896 matrices of size 4n-by-4n by first sending each quaternion entry
897 `z = a + bi + cj + dk` to the block-complex matrix
898 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
903 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
907 sage: Q = QuaternionAlgebra(QQ,-1,-1)
908 sage: i,j,k = Q.gens()
909 sage: x = 1 + 2*i + 3*j + 4*k
910 sage: M = matrix(Q, 1, [[x]])
911 sage: _embed_quaternion_matrix(M)
917 Embedding is a homomorphism (isomorphism, in fact)::
919 sage: set_random_seed()
920 sage: n = ZZ.random_element(5)
921 sage: Q = QuaternionAlgebra(QQ,-1,-1)
922 sage: X = random_matrix(Q, n)
923 sage: Y = random_matrix(Q, n)
924 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
925 sage: expected = _embed_quaternion_matrix(X*Y)
926 sage: actual == expected
930 quaternions
= M
.base_ring()
933 raise ValueError("the matrix 'M' must be square")
935 F
= QuadraticField(-1, 'i')
940 t
= z
.coefficient_tuple()
945 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
946 [-c
+ d
*i
, a
- b
*i
]])
947 blocks
.append(_embed_complex_matrix(cplx_matrix
))
949 # We should have real entries by now, so use the realest field
950 # we've got for the return value.
951 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
954 def _unembed_quaternion_matrix(M
):
956 The inverse of _embed_quaternion_matrix().
960 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
961 ....: _unembed_quaternion_matrix)
965 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
966 ....: [-2, 1, -4, 3],
967 ....: [-3, 4, 1, -2],
968 ....: [-4, -3, 2, 1]])
969 sage: _unembed_quaternion_matrix(M)
970 [1 + 2*i + 3*j + 4*k]
974 Unembedding is the inverse of embedding::
976 sage: set_random_seed()
977 sage: Q = QuaternionAlgebra(QQ, -1, -1)
978 sage: M = random_matrix(Q, 3)
979 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
985 raise ValueError("the matrix 'M' must be square")
986 if not n
.mod(4).is_zero():
987 raise ValueError("the matrix 'M' must be a complex embedding")
989 Q
= QuaternionAlgebra(QQ
,-1,-1)
992 # Go top-left to bottom-right (reading order), converting every
993 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
996 for l
in xrange(n
/4):
997 for m
in xrange(n
/4):
998 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
999 if submat
[0,0] != submat
[1,1].conjugate():
1000 raise ValueError('bad on-diagonal submatrix')
1001 if submat
[0,1] != -submat
[1,0].conjugate():
1002 raise ValueError('bad off-diagonal submatrix')
1003 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1004 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1007 return matrix(Q
, n
/4, elements
)
1010 # The usual inner product on R^n.
1012 return x
.vector().inner_product(y
.vector())
1014 # The inner product used for the real symmetric simple EJA.
1015 # We keep it as a separate function because e.g. the complex
1016 # algebra uses the same inner product, except divided by 2.
1017 def _matrix_ip(X
,Y
):
1018 X_mat
= X
.natural_representation()
1019 Y_mat
= Y
.natural_representation()
1020 return (X_mat
*Y_mat
).trace()
1023 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1025 The rank-n simple EJA consisting of real symmetric n-by-n
1026 matrices, the usual symmetric Jordan product, and the trace inner
1027 product. It has dimension `(n^2 + n)/2` over the reals.
1031 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1035 sage: J = RealSymmetricEJA(2)
1036 sage: e0, e1, e2 = J.gens()
1046 The degree of this algebra is `(n^2 + n) / 2`::
1048 sage: set_random_seed()
1049 sage: n = ZZ.random_element(1,5)
1050 sage: J = RealSymmetricEJA(n)
1051 sage: J.degree() == (n^2 + n)/2
1054 The Jordan multiplication is what we think it is::
1056 sage: set_random_seed()
1057 sage: n = ZZ.random_element(1,5)
1058 sage: J = RealSymmetricEJA(n)
1059 sage: x = J.random_element()
1060 sage: y = J.random_element()
1061 sage: actual = (x*y).natural_representation()
1062 sage: X = x.natural_representation()
1063 sage: Y = y.natural_representation()
1064 sage: expected = (X*Y + Y*X)/2
1065 sage: actual == expected
1067 sage: J(expected) == x*y
1072 def __classcall_private__(cls
, n
, field
=QQ
):
1073 S
= _real_symmetric_basis(n
, field
=field
)
1074 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
1076 fdeja
= super(RealSymmetricEJA
, cls
)
1077 return fdeja
.__classcall
_private
__(cls
,
1083 def inner_product(self
, x
, y
):
1084 return _matrix_ip(x
,y
)
1087 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1089 The rank-n simple EJA consisting of complex Hermitian n-by-n
1090 matrices over the real numbers, the usual symmetric Jordan product,
1091 and the real-part-of-trace inner product. It has dimension `n^2` over
1096 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1100 The degree of this algebra is `n^2`::
1102 sage: set_random_seed()
1103 sage: n = ZZ.random_element(1,5)
1104 sage: J = ComplexHermitianEJA(n)
1105 sage: J.degree() == n^2
1108 The Jordan multiplication is what we think it is::
1110 sage: set_random_seed()
1111 sage: n = ZZ.random_element(1,5)
1112 sage: J = ComplexHermitianEJA(n)
1113 sage: x = J.random_element()
1114 sage: y = J.random_element()
1115 sage: actual = (x*y).natural_representation()
1116 sage: X = x.natural_representation()
1117 sage: Y = y.natural_representation()
1118 sage: expected = (X*Y + Y*X)/2
1119 sage: actual == expected
1121 sage: J(expected) == x*y
1126 def __classcall_private__(cls
, n
, field
=QQ
):
1127 S
= _complex_hermitian_basis(n
)
1128 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
1130 fdeja
= super(ComplexHermitianEJA
, cls
)
1131 return fdeja
.__classcall
_private
__(cls
,
1137 def inner_product(self
, x
, y
):
1138 # Since a+bi on the diagonal is represented as
1143 # we'll double-count the "a" entries if we take the trace of
1145 return _matrix_ip(x
,y
)/2
1148 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1150 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1151 matrices, the usual symmetric Jordan product, and the
1152 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1157 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1161 The degree of this algebra is `n^2`::
1163 sage: set_random_seed()
1164 sage: n = ZZ.random_element(1,5)
1165 sage: J = QuaternionHermitianEJA(n)
1166 sage: J.degree() == 2*(n^2) - n
1169 The Jordan multiplication is what we think it is::
1171 sage: set_random_seed()
1172 sage: n = ZZ.random_element(1,5)
1173 sage: J = QuaternionHermitianEJA(n)
1174 sage: x = J.random_element()
1175 sage: y = J.random_element()
1176 sage: actual = (x*y).natural_representation()
1177 sage: X = x.natural_representation()
1178 sage: Y = y.natural_representation()
1179 sage: expected = (X*Y + Y*X)/2
1180 sage: actual == expected
1182 sage: J(expected) == x*y
1187 def __classcall_private__(cls
, n
, field
=QQ
):
1188 S
= _quaternion_hermitian_basis(n
)
1189 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
1191 fdeja
= super(QuaternionHermitianEJA
, cls
)
1192 return fdeja
.__classcall
_private
__(cls
,
1198 def inner_product(self
, x
, y
):
1199 # Since a+bi+cj+dk on the diagonal is represented as
1201 # a + bi +cj + dk = [ a b c d]
1206 # we'll quadruple-count the "a" entries if we take the trace of
1208 return _matrix_ip(x
,y
)/4
1211 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1213 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1214 with the usual inner product and jordan product ``x*y =
1215 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1220 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1224 This multiplication table can be verified by hand::
1226 sage: J = JordanSpinEJA(4)
1227 sage: e0,e1,e2,e3 = J.gens()
1245 def __classcall_private__(cls
, n
, field
=QQ
):
1247 id_matrix
= matrix
.identity(field
, n
)
1249 ei
= id_matrix
.column(i
)
1250 Qi
= matrix
.zero(field
, n
)
1252 Qi
.set_column(0, ei
)
1253 Qi
+= matrix
.diagonal(n
, [ei
[0]]*n
)
1254 # The addition of the diagonal matrix adds an extra ei[0] in the
1255 # upper-left corner of the matrix.
1256 Qi
[0,0] = Qi
[0,0] * ~
field(2)
1259 # The rank of the spin algebra is two, unless we're in a
1260 # one-dimensional ambient space (because the rank is bounded by
1261 # the ambient dimension).
1262 fdeja
= super(JordanSpinEJA
, cls
)
1263 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=min(n
,2))
1265 def inner_product(self
, x
, y
):
1266 return _usual_ip(x
,y
)