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.
235 l1
= [matrix
.column(W
.coordinates((x
**k
).vector())) for k
in range(r
)]
236 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
237 A_of_x
= matrix
.block(R
, 1, n
, (l1
+ l2
))
238 xr
= W
.coordinates((x
**r
).vector())
239 return (A_of_x
, x
, xr
, A_of_x
.det())
243 def characteristic_polynomial(self
):
245 Return a characteristic polynomial that works for all elements
248 The resulting polynomial has `n+1` variables, where `n` is the
249 dimension of this algebra. The first `n` variables correspond to
250 the coordinates of an algebra element: when evaluated at the
251 coordinates of an algebra element with respect to a certain
252 basis, the result is a univariate polynomial (in the one
253 remaining variable ``t``), namely the characteristic polynomial
258 sage: from mjo.eja.eja_algebra import JordanSpinEJA
262 The characteristic polynomial in the spin algebra is given in
263 Alizadeh, Example 11.11::
265 sage: J = JordanSpinEJA(3)
266 sage: p = J.characteristic_polynomial(); p
267 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
268 sage: xvec = J.one().vector()
276 # The list of coefficient polynomials a_1, a_2, ..., a_n.
277 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
279 # We go to a bit of trouble here to reorder the
280 # indeterminates, so that it's easier to evaluate the
281 # characteristic polynomial at x's coordinates and get back
282 # something in terms of t, which is what we want.
284 S
= PolynomialRing(self
.base_ring(),'t')
286 S
= PolynomialRing(S
, R
.variable_names())
289 # Note: all entries past the rth should be zero. The
290 # coefficient of the highest power (x^r) is 1, but it doesn't
291 # appear in the solution vector which contains coefficients
292 # for the other powers (to make them sum to x^r).
294 a
[r
] = 1 # corresponds to x^r
296 # When the rank is equal to the dimension, trying to
297 # assign a[r] goes out-of-bounds.
298 a
.append(1) # corresponds to x^r
300 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
303 def inner_product(self
, x
, y
):
305 The inner product associated with this Euclidean Jordan algebra.
307 Defaults to the trace inner product, but can be overridden by
308 subclasses if they are sure that the necessary properties are
313 sage: from mjo.eja.eja_algebra import random_eja
317 The inner product must satisfy its axiom for this algebra to truly
318 be a Euclidean Jordan Algebra::
320 sage: set_random_seed()
321 sage: J = random_eja()
322 sage: x = J.random_element()
323 sage: y = J.random_element()
324 sage: z = J.random_element()
325 sage: (x*y).inner_product(z) == y.inner_product(x*z)
329 if (not x
in self
) or (not y
in self
):
330 raise TypeError("arguments must live in this algebra")
331 return x
.trace_inner_product(y
)
334 def natural_basis(self
):
336 Return a more-natural representation of this algebra's basis.
338 Every finite-dimensional Euclidean Jordan Algebra is a direct
339 sum of five simple algebras, four of which comprise Hermitian
340 matrices. This method returns the original "natural" basis
341 for our underlying vector space. (Typically, the natural basis
342 is used to construct the multiplication table in the first place.)
344 Note that this will always return a matrix. The standard basis
345 in `R^n` will be returned as `n`-by-`1` column matrices.
349 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
350 ....: RealSymmetricEJA)
354 sage: J = RealSymmetricEJA(2)
357 sage: J.natural_basis()
365 sage: J = JordanSpinEJA(2)
368 sage: J.natural_basis()
375 if self
._natural
_basis
is None:
376 return tuple( b
.vector().column() for b
in self
.basis() )
378 return self
._natural
_basis
383 Return the rank of this EJA.
387 The author knows of no algorithm to compute the rank of an EJA
388 where only the multiplication table is known. In lieu of one, we
389 require the rank to be specified when the algebra is created,
390 and simply pass along that number here.
394 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
395 ....: RealSymmetricEJA,
396 ....: ComplexHermitianEJA,
397 ....: QuaternionHermitianEJA,
402 The rank of the Jordan spin algebra is always two::
404 sage: JordanSpinEJA(2).rank()
406 sage: JordanSpinEJA(3).rank()
408 sage: JordanSpinEJA(4).rank()
411 The rank of the `n`-by-`n` Hermitian real, complex, or
412 quaternion matrices is `n`::
414 sage: RealSymmetricEJA(2).rank()
416 sage: ComplexHermitianEJA(2).rank()
418 sage: QuaternionHermitianEJA(2).rank()
420 sage: RealSymmetricEJA(5).rank()
422 sage: ComplexHermitianEJA(5).rank()
424 sage: QuaternionHermitianEJA(5).rank()
429 Ensure that every EJA that we know how to construct has a
430 positive integer rank::
432 sage: set_random_seed()
433 sage: r = random_eja().rank()
434 sage: r in ZZ and r > 0
441 def vector_space(self
):
443 Return the vector space that underlies this algebra.
447 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
451 sage: J = RealSymmetricEJA(2)
452 sage: J.vector_space()
453 Vector space of dimension 3 over Rational Field
456 return self
.zero().vector().parent().ambient_vector_space()
459 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
462 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
464 Return the Euclidean Jordan Algebra corresponding to the set
465 `R^n` under the Hadamard product.
467 Note: this is nothing more than the Cartesian product of ``n``
468 copies of the spin algebra. Once Cartesian product algebras
469 are implemented, this can go.
473 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
477 This multiplication table can be verified by hand::
479 sage: J = RealCartesianProductEJA(3)
480 sage: e0,e1,e2 = J.gens()
496 def __classcall_private__(cls
, n
, field
=QQ
):
497 # The FiniteDimensionalAlgebra constructor takes a list of
498 # matrices, the ith representing right multiplication by the ith
499 # basis element in the vector space. So if e_1 = (1,0,0), then
500 # right (Hadamard) multiplication of x by e_1 picks out the first
501 # component of x; and likewise for the ith basis element e_i.
502 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
505 fdeja
= super(RealCartesianProductEJA
, cls
)
506 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=n
)
508 def inner_product(self
, x
, y
):
509 return _usual_ip(x
,y
)
514 Return a "random" finite-dimensional Euclidean Jordan Algebra.
518 For now, we choose a random natural number ``n`` (greater than zero)
519 and then give you back one of the following:
521 * The cartesian product of the rational numbers ``n`` times; this is
522 ``QQ^n`` with the Hadamard product.
524 * The Jordan spin algebra on ``QQ^n``.
526 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
529 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
530 in the space of ``2n``-by-``2n`` real symmetric matrices.
532 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
533 in the space of ``4n``-by-``4n`` real symmetric matrices.
535 Later this might be extended to return Cartesian products of the
540 sage: from mjo.eja.eja_algebra import random_eja
545 Euclidean Jordan algebra of degree...
549 # The max_n component lets us choose different upper bounds on the
550 # value "n" that gets passed to the constructor. This is needed
551 # because e.g. R^{10} is reasonable to test, while the Hermitian
552 # 10-by-10 quaternion matrices are not.
553 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
555 (RealSymmetricEJA
, 5),
556 (ComplexHermitianEJA
, 4),
557 (QuaternionHermitianEJA
, 3)])
558 n
= ZZ
.random_element(1, max_n
)
559 return constructor(n
, field
=QQ
)
563 def _real_symmetric_basis(n
, field
=QQ
):
565 Return a basis for the space of real symmetric n-by-n matrices.
567 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
571 for j
in xrange(i
+1):
572 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
576 # Beware, orthogonal but not normalized!
577 Sij
= Eij
+ Eij
.transpose()
582 def _complex_hermitian_basis(n
, field
=QQ
):
584 Returns a basis for the space of complex Hermitian n-by-n matrices.
588 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
592 sage: set_random_seed()
593 sage: n = ZZ.random_element(1,5)
594 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
598 F
= QuadraticField(-1, 'I')
601 # This is like the symmetric case, but we need to be careful:
603 # * We want conjugate-symmetry, not just symmetry.
604 # * The diagonal will (as a result) be real.
608 for j
in xrange(i
+1):
609 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
611 Sij
= _embed_complex_matrix(Eij
)
614 # Beware, orthogonal but not normalized! The second one
615 # has a minus because it's conjugated.
616 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
618 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
623 def _quaternion_hermitian_basis(n
, field
=QQ
):
625 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
629 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
633 sage: set_random_seed()
634 sage: n = ZZ.random_element(1,5)
635 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
639 Q
= QuaternionAlgebra(QQ
,-1,-1)
642 # This is like the symmetric case, but we need to be careful:
644 # * We want conjugate-symmetry, not just symmetry.
645 # * The diagonal will (as a result) be real.
649 for j
in xrange(i
+1):
650 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
652 Sij
= _embed_quaternion_matrix(Eij
)
655 # Beware, orthogonal but not normalized! The second,
656 # third, and fourth ones have a minus because they're
658 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
660 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
662 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
664 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
670 def _multiplication_table_from_matrix_basis(basis
):
672 At least three of the five simple Euclidean Jordan algebras have the
673 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
674 multiplication on the right is matrix multiplication. Given a basis
675 for the underlying matrix space, this function returns a
676 multiplication table (obtained by looping through the basis
677 elements) for an algebra of those matrices. A reordered copy
678 of the basis is also returned to work around the fact that
679 the ``span()`` in this function will change the order of the basis
680 from what we think it is, to... something else.
682 # In S^2, for example, we nominally have four coordinates even
683 # though the space is of dimension three only. The vector space V
684 # is supposed to hold the entire long vector, and the subspace W
685 # of V will be spanned by the vectors that arise from symmetric
686 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
687 field
= basis
[0].base_ring()
688 dimension
= basis
[0].nrows()
690 V
= VectorSpace(field
, dimension
**2)
691 W
= V
.span( _mat2vec(s
) for s
in basis
)
693 # Taking the span above reorders our basis (thanks, jerk!) so we
694 # need to put our "matrix basis" in the same order as the
695 # (reordered) vector basis.
696 S
= tuple( _vec2mat(b
) for b
in W
.basis() )
700 # Brute force the multiplication-by-s matrix by looping
701 # through all elements of the basis and doing the computation
702 # to find out what the corresponding row should be. BEWARE:
703 # these multiplication tables won't be symmetric! It therefore
704 # becomes REALLY IMPORTANT that the underlying algebra
705 # constructor uses ROW vectors and not COLUMN vectors. That's
706 # why we're computing rows here and not columns.
709 this_row
= _mat2vec((s
*t
+ t
*s
)/2)
710 Q_rows
.append(W
.coordinates(this_row
))
711 Q
= matrix(field
, W
.dimension(), Q_rows
)
717 def _embed_complex_matrix(M
):
719 Embed the n-by-n complex matrix ``M`` into the space of real
720 matrices of size 2n-by-2n via the map the sends each entry `z = a +
721 bi` to the block matrix ``[[a,b],[-b,a]]``.
725 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
729 sage: F = QuadraticField(-1,'i')
730 sage: x1 = F(4 - 2*i)
731 sage: x2 = F(1 + 2*i)
734 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
735 sage: _embed_complex_matrix(M)
744 Embedding is a homomorphism (isomorphism, in fact)::
746 sage: set_random_seed()
747 sage: n = ZZ.random_element(5)
748 sage: F = QuadraticField(-1, 'i')
749 sage: X = random_matrix(F, n)
750 sage: Y = random_matrix(F, n)
751 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
752 sage: expected = _embed_complex_matrix(X*Y)
753 sage: actual == expected
759 raise ValueError("the matrix 'M' must be square")
760 field
= M
.base_ring()
765 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
767 # We can drop the imaginaries here.
768 return matrix
.block(field
.base_ring(), n
, blocks
)
771 def _unembed_complex_matrix(M
):
773 The inverse of _embed_complex_matrix().
777 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
778 ....: _unembed_complex_matrix)
782 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
783 ....: [-2, 1, -4, 3],
784 ....: [ 9, 10, 11, 12],
785 ....: [-10, 9, -12, 11] ])
786 sage: _unembed_complex_matrix(A)
788 [ 10*i + 9 12*i + 11]
792 Unembedding is the inverse of embedding::
794 sage: set_random_seed()
795 sage: F = QuadraticField(-1, 'i')
796 sage: M = random_matrix(F, 3)
797 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
803 raise ValueError("the matrix 'M' must be square")
804 if not n
.mod(2).is_zero():
805 raise ValueError("the matrix 'M' must be a complex embedding")
807 F
= QuadraticField(-1, 'i')
810 # Go top-left to bottom-right (reading order), converting every
811 # 2-by-2 block we see to a single complex element.
813 for k
in xrange(n
/2):
814 for j
in xrange(n
/2):
815 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
816 if submat
[0,0] != submat
[1,1]:
817 raise ValueError('bad on-diagonal submatrix')
818 if submat
[0,1] != -submat
[1,0]:
819 raise ValueError('bad off-diagonal submatrix')
820 z
= submat
[0,0] + submat
[0,1]*i
823 return matrix(F
, n
/2, elements
)
826 def _embed_quaternion_matrix(M
):
828 Embed the n-by-n quaternion matrix ``M`` into the space of real
829 matrices of size 4n-by-4n by first sending each quaternion entry
830 `z = a + bi + cj + dk` to the block-complex matrix
831 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
836 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
840 sage: Q = QuaternionAlgebra(QQ,-1,-1)
841 sage: i,j,k = Q.gens()
842 sage: x = 1 + 2*i + 3*j + 4*k
843 sage: M = matrix(Q, 1, [[x]])
844 sage: _embed_quaternion_matrix(M)
850 Embedding is a homomorphism (isomorphism, in fact)::
852 sage: set_random_seed()
853 sage: n = ZZ.random_element(5)
854 sage: Q = QuaternionAlgebra(QQ,-1,-1)
855 sage: X = random_matrix(Q, n)
856 sage: Y = random_matrix(Q, n)
857 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
858 sage: expected = _embed_quaternion_matrix(X*Y)
859 sage: actual == expected
863 quaternions
= M
.base_ring()
866 raise ValueError("the matrix 'M' must be square")
868 F
= QuadraticField(-1, 'i')
873 t
= z
.coefficient_tuple()
878 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
879 [-c
+ d
*i
, a
- b
*i
]])
880 blocks
.append(_embed_complex_matrix(cplx_matrix
))
882 # We should have real entries by now, so use the realest field
883 # we've got for the return value.
884 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
887 def _unembed_quaternion_matrix(M
):
889 The inverse of _embed_quaternion_matrix().
893 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
894 ....: _unembed_quaternion_matrix)
898 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
899 ....: [-2, 1, -4, 3],
900 ....: [-3, 4, 1, -2],
901 ....: [-4, -3, 2, 1]])
902 sage: _unembed_quaternion_matrix(M)
903 [1 + 2*i + 3*j + 4*k]
907 Unembedding is the inverse of embedding::
909 sage: set_random_seed()
910 sage: Q = QuaternionAlgebra(QQ, -1, -1)
911 sage: M = random_matrix(Q, 3)
912 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
918 raise ValueError("the matrix 'M' must be square")
919 if not n
.mod(4).is_zero():
920 raise ValueError("the matrix 'M' must be a complex embedding")
922 Q
= QuaternionAlgebra(QQ
,-1,-1)
925 # Go top-left to bottom-right (reading order), converting every
926 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
929 for l
in xrange(n
/4):
930 for m
in xrange(n
/4):
931 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
932 if submat
[0,0] != submat
[1,1].conjugate():
933 raise ValueError('bad on-diagonal submatrix')
934 if submat
[0,1] != -submat
[1,0].conjugate():
935 raise ValueError('bad off-diagonal submatrix')
936 z
= submat
[0,0].real() + submat
[0,0].imag()*i
937 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
940 return matrix(Q
, n
/4, elements
)
943 # The usual inner product on R^n.
945 return x
.vector().inner_product(y
.vector())
947 # The inner product used for the real symmetric simple EJA.
948 # We keep it as a separate function because e.g. the complex
949 # algebra uses the same inner product, except divided by 2.
951 X_mat
= X
.natural_representation()
952 Y_mat
= Y
.natural_representation()
953 return (X_mat
*Y_mat
).trace()
956 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
958 The rank-n simple EJA consisting of real symmetric n-by-n
959 matrices, the usual symmetric Jordan product, and the trace inner
960 product. It has dimension `(n^2 + n)/2` over the reals.
964 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
968 sage: J = RealSymmetricEJA(2)
969 sage: e0, e1, e2 = J.gens()
979 The degree of this algebra is `(n^2 + n) / 2`::
981 sage: set_random_seed()
982 sage: n = ZZ.random_element(1,5)
983 sage: J = RealSymmetricEJA(n)
984 sage: J.degree() == (n^2 + n)/2
987 The Jordan multiplication is what we think it is::
989 sage: set_random_seed()
990 sage: n = ZZ.random_element(1,5)
991 sage: J = RealSymmetricEJA(n)
992 sage: x = J.random_element()
993 sage: y = J.random_element()
994 sage: actual = (x*y).natural_representation()
995 sage: X = x.natural_representation()
996 sage: Y = y.natural_representation()
997 sage: expected = (X*Y + Y*X)/2
998 sage: actual == expected
1000 sage: J(expected) == x*y
1005 def __classcall_private__(cls
, n
, field
=QQ
):
1006 S
= _real_symmetric_basis(n
, field
=field
)
1007 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
1009 fdeja
= super(RealSymmetricEJA
, cls
)
1010 return fdeja
.__classcall
_private
__(cls
,
1016 def inner_product(self
, x
, y
):
1017 return _matrix_ip(x
,y
)
1020 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1022 The rank-n simple EJA consisting of complex Hermitian n-by-n
1023 matrices over the real numbers, the usual symmetric Jordan product,
1024 and the real-part-of-trace inner product. It has dimension `n^2` over
1029 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1033 The degree of this algebra is `n^2`::
1035 sage: set_random_seed()
1036 sage: n = ZZ.random_element(1,5)
1037 sage: J = ComplexHermitianEJA(n)
1038 sage: J.degree() == n^2
1041 The Jordan multiplication is what we think it is::
1043 sage: set_random_seed()
1044 sage: n = ZZ.random_element(1,5)
1045 sage: J = ComplexHermitianEJA(n)
1046 sage: x = J.random_element()
1047 sage: y = J.random_element()
1048 sage: actual = (x*y).natural_representation()
1049 sage: X = x.natural_representation()
1050 sage: Y = y.natural_representation()
1051 sage: expected = (X*Y + Y*X)/2
1052 sage: actual == expected
1054 sage: J(expected) == x*y
1059 def __classcall_private__(cls
, n
, field
=QQ
):
1060 S
= _complex_hermitian_basis(n
)
1061 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
1063 fdeja
= super(ComplexHermitianEJA
, cls
)
1064 return fdeja
.__classcall
_private
__(cls
,
1070 def inner_product(self
, x
, y
):
1071 # Since a+bi on the diagonal is represented as
1076 # we'll double-count the "a" entries if we take the trace of
1078 return _matrix_ip(x
,y
)/2
1081 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1083 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1084 matrices, the usual symmetric Jordan product, and the
1085 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1090 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1094 The degree of this algebra is `n^2`::
1096 sage: set_random_seed()
1097 sage: n = ZZ.random_element(1,5)
1098 sage: J = QuaternionHermitianEJA(n)
1099 sage: J.degree() == 2*(n^2) - n
1102 The Jordan multiplication is what we think it is::
1104 sage: set_random_seed()
1105 sage: n = ZZ.random_element(1,5)
1106 sage: J = QuaternionHermitianEJA(n)
1107 sage: x = J.random_element()
1108 sage: y = J.random_element()
1109 sage: actual = (x*y).natural_representation()
1110 sage: X = x.natural_representation()
1111 sage: Y = y.natural_representation()
1112 sage: expected = (X*Y + Y*X)/2
1113 sage: actual == expected
1115 sage: J(expected) == x*y
1120 def __classcall_private__(cls
, n
, field
=QQ
):
1121 S
= _quaternion_hermitian_basis(n
)
1122 (Qs
, T
) = _multiplication_table_from_matrix_basis(S
)
1124 fdeja
= super(QuaternionHermitianEJA
, cls
)
1125 return fdeja
.__classcall
_private
__(cls
,
1131 def inner_product(self
, x
, y
):
1132 # Since a+bi+cj+dk on the diagonal is represented as
1134 # a + bi +cj + dk = [ a b c d]
1139 # we'll quadruple-count the "a" entries if we take the trace of
1141 return _matrix_ip(x
,y
)/4
1144 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1146 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1147 with the usual inner product and jordan product ``x*y =
1148 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1153 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1157 This multiplication table can be verified by hand::
1159 sage: J = JordanSpinEJA(4)
1160 sage: e0,e1,e2,e3 = J.gens()
1178 def __classcall_private__(cls
, n
, field
=QQ
):
1180 id_matrix
= matrix
.identity(field
, n
)
1182 ei
= id_matrix
.column(i
)
1183 Qi
= matrix
.zero(field
, n
)
1185 Qi
.set_column(0, ei
)
1186 Qi
+= matrix
.diagonal(n
, [ei
[0]]*n
)
1187 # The addition of the diagonal matrix adds an extra ei[0] in the
1188 # upper-left corner of the matrix.
1189 Qi
[0,0] = Qi
[0,0] * ~
field(2)
1192 # The rank of the spin algebra is two, unless we're in a
1193 # one-dimensional ambient space (because the rank is bounded by
1194 # the ambient dimension).
1195 fdeja
= super(JordanSpinEJA
, cls
)
1196 return fdeja
.__classcall
_private
__(cls
, field
, Qs
, rank
=min(n
,2))
1198 def inner_product(self
, x
, y
):
1199 return _usual_ip(x
,y
)