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
.quatalg
.quaternion_algebra
import QuaternionAlgebra
9 from sage
.categories
.finite_dimensional_algebras_with_basis
import FiniteDimensionalAlgebrasWithBasis
10 from sage
.combinat
.free_module
import CombinatorialFreeModule
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
21 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
22 from mjo
.eja
.eja_utils
import _mat2vec
24 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
35 sage: from mjo.eja.eja_algebra import random_eja
39 By definition, Jordan multiplication commutes::
41 sage: set_random_seed()
42 sage: J = random_eja()
43 sage: x = J.random_element()
44 sage: y = J.random_element()
50 self
._natural
_basis
= natural_basis
53 category
= FiniteDimensionalAlgebrasWithBasis(field
).Unital()
54 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
56 range(len(mult_table
)),
59 self
.print_options(bracket
='')
61 # The multiplication table we're given is necessarily in terms
62 # of vectors, because we don't have an algebra yet for
63 # anything to be an element of. However, it's faster in the
64 # long run to have the multiplication table be in terms of
65 # algebra elements. We do this after calling the superclass
66 # constructor so that from_vector() knows what to do.
67 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
68 for ls
in mult_table
]
71 def _element_constructor_(self
, elt
):
73 Construct an element of this algebra from its natural
76 This gets called only after the parent element _call_ method
77 fails to find a coercion for the argument.
81 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
82 ....: RealCartesianProductEJA,
83 ....: RealSymmetricEJA)
87 The identity in `S^n` is converted to the identity in the EJA::
89 sage: J = RealSymmetricEJA(3)
90 sage: I = matrix.identity(QQ,3)
94 This skew-symmetric matrix can't be represented in the EJA::
96 sage: J = RealSymmetricEJA(3)
97 sage: A = matrix(QQ,3, lambda i,j: i-j)
99 Traceback (most recent call last):
101 ArithmeticError: vector is not in free module
105 Ensure that we can convert any element of the two non-matrix
106 simple algebras (whose natural representations are their usual
107 vector representations) back and forth faithfully::
109 sage: set_random_seed()
110 sage: J = RealCartesianProductEJA(5)
111 sage: x = J.random_element()
112 sage: J(x.to_vector().column()) == x
114 sage: J = JordanSpinEJA(5)
115 sage: x = J.random_element()
116 sage: J(x.to_vector().column()) == x
120 natural_basis
= self
.natural_basis()
121 if elt
not in natural_basis
[0].matrix_space():
122 raise ValueError("not a naturally-represented algebra element")
124 # Thanks for nothing! Matrix spaces aren't vector
125 # spaces in Sage, so we have to figure out its
126 # natural-basis coordinates ourselves.
127 V
= VectorSpace(elt
.base_ring(), elt
.nrows()*elt
.ncols())
128 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
129 coords
= W
.coordinate_vector(_mat2vec(elt
))
130 return self
.from_vector(coords
)
135 Return a string representation of ``self``.
139 sage: from mjo.eja.eja_algebra import JordanSpinEJA
143 Ensure that it says what we think it says::
145 sage: JordanSpinEJA(2, field=QQ)
146 Euclidean Jordan algebra of degree 2 over Rational Field
147 sage: JordanSpinEJA(3, field=RDF)
148 Euclidean Jordan algebra of degree 3 over Real Double Field
151 # TODO: change this to say "dimension" and fix all the tests.
152 fmt
= "Euclidean Jordan algebra of degree {} over {}"
153 return fmt
.format(self
.dimension(), self
.base_ring())
155 def product_on_basis(self
, i
, j
):
156 return self
._multiplication
_table
[i
][j
]
158 def _a_regular_element(self
):
160 Guess a regular element. Needed to compute the basis for our
161 characteristic polynomial coefficients.
165 sage: from mjo.eja.eja_algebra import random_eja
169 Ensure that this hacky method succeeds for every algebra that we
170 know how to construct::
172 sage: set_random_seed()
173 sage: J = random_eja()
174 sage: J._a_regular_element().is_regular()
179 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
180 if not z
.is_regular():
181 raise ValueError("don't know a regular element")
186 def _charpoly_basis_space(self
):
188 Return the vector space spanned by the basis used in our
189 characteristic polynomial coefficients. This is used not only to
190 compute those coefficients, but also any time we need to
191 evaluate the coefficients (like when we compute the trace or
194 z
= self
._a
_regular
_element
()
195 V
= self
.vector_space()
196 V1
= V
.span_of_basis( (z
**k
).to_vector() for k
in range(self
.rank()) )
197 b
= (V1
.basis() + V1
.complement().basis())
198 return V
.span_of_basis(b
)
202 def _charpoly_coeff(self
, i
):
204 Return the coefficient polynomial "a_{i}" of this algebra's
205 general characteristic polynomial.
207 Having this be a separate cached method lets us compute and
208 store the trace/determinant (a_{r-1} and a_{0} respectively)
209 separate from the entire characteristic polynomial.
211 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
212 R
= A_of_x
.base_ring()
214 # Guaranteed by theory
217 # Danger: the in-place modification is done for performance
218 # reasons (reconstructing a matrix with huge polynomial
219 # entries is slow), but I don't know how cached_method works,
220 # so it's highly possible that we're modifying some global
221 # list variable by reference, here. In other words, you
222 # probably shouldn't call this method twice on the same
223 # algebra, at the same time, in two threads
224 Ai_orig
= A_of_x
.column(i
)
225 A_of_x
.set_column(i
,xr
)
226 numerator
= A_of_x
.det()
227 A_of_x
.set_column(i
,Ai_orig
)
229 # We're relying on the theory here to ensure that each a_i is
230 # indeed back in R, and the added negative signs are to make
231 # the whole charpoly expression sum to zero.
232 return R(-numerator
/detA
)
236 def _charpoly_matrix_system(self
):
238 Compute the matrix whose entries A_ij are polynomials in
239 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
240 corresponding to `x^r` and the determinent of the matrix A =
241 [A_ij]. In other words, all of the fixed (cachable) data needed
242 to compute the coefficients of the characteristic polynomial.
247 # Construct a new algebra over a multivariate polynomial ring...
248 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
249 R
= PolynomialRing(self
.base_ring(), names
)
250 # Hack around the fact that our multiplication table is in terms of
251 # algebra elements but the constructor wants it in terms of vectors.
252 vmt
= [ tuple(map(lambda x
: x
.to_vector(), ls
))
253 for ls
in self
._multiplication
_table
]
254 J
= FiniteDimensionalEuclideanJordanAlgebra(R
, tuple(vmt
), r
)
256 idmat
= matrix
.identity(J
.base_ring(), n
)
258 W
= self
._charpoly
_basis
_space
()
259 W
= W
.change_ring(R
.fraction_field())
261 # Starting with the standard coordinates x = (X1,X2,...,Xn)
262 # and then converting the entries to W-coordinates allows us
263 # to pass in the standard coordinates to the charpoly and get
264 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
267 # W.coordinates(x^2) eval'd at (standard z-coords)
271 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
273 # We want the middle equivalent thing in our matrix, but use
274 # the first equivalent thing instead so that we can pass in
275 # standard coordinates.
276 x
= J
.from_vector(W(R
.gens()))
278 # Handle the zeroth power separately, because computing
279 # the unit element in J is mathematically suspect.
280 x0
= W
.coordinate_vector(self
.one().to_vector())
282 l1
+= [ W
.coordinate_vector((x
**k
).to_vector()).column()
283 for k
in range(1,r
) ]
284 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
285 A_of_x
= matrix
.block(R
, 1, n
, (l1
+ l2
))
286 xr
= W
.coordinate_vector((x
**r
).to_vector())
287 return (A_of_x
, x
, xr
, A_of_x
.det())
291 def characteristic_polynomial(self
):
293 Return a characteristic polynomial that works for all elements
296 The resulting polynomial has `n+1` variables, where `n` is the
297 dimension of this algebra. The first `n` variables correspond to
298 the coordinates of an algebra element: when evaluated at the
299 coordinates of an algebra element with respect to a certain
300 basis, the result is a univariate polynomial (in the one
301 remaining variable ``t``), namely the characteristic polynomial
306 sage: from mjo.eja.eja_algebra import JordanSpinEJA
310 The characteristic polynomial in the spin algebra is given in
311 Alizadeh, Example 11.11::
313 sage: J = JordanSpinEJA(3)
314 sage: p = J.characteristic_polynomial(); p
315 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
316 sage: xvec = J.one().to_vector()
324 # The list of coefficient polynomials a_1, a_2, ..., a_n.
325 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
327 # We go to a bit of trouble here to reorder the
328 # indeterminates, so that it's easier to evaluate the
329 # characteristic polynomial at x's coordinates and get back
330 # something in terms of t, which is what we want.
332 S
= PolynomialRing(self
.base_ring(),'t')
334 S
= PolynomialRing(S
, R
.variable_names())
337 # Note: all entries past the rth should be zero. The
338 # coefficient of the highest power (x^r) is 1, but it doesn't
339 # appear in the solution vector which contains coefficients
340 # for the other powers (to make them sum to x^r).
342 a
[r
] = 1 # corresponds to x^r
344 # When the rank is equal to the dimension, trying to
345 # assign a[r] goes out-of-bounds.
346 a
.append(1) # corresponds to x^r
348 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
351 def inner_product(self
, x
, y
):
353 The inner product associated with this Euclidean Jordan algebra.
355 Defaults to the trace inner product, but can be overridden by
356 subclasses if they are sure that the necessary properties are
361 sage: from mjo.eja.eja_algebra import random_eja
365 The inner product must satisfy its axiom for this algebra to truly
366 be a Euclidean Jordan Algebra::
368 sage: set_random_seed()
369 sage: J = random_eja()
370 sage: x = J.random_element()
371 sage: y = J.random_element()
372 sage: z = J.random_element()
373 sage: (x*y).inner_product(z) == y.inner_product(x*z)
377 if (not x
in self
) or (not y
in self
):
378 raise TypeError("arguments must live in this algebra")
379 return x
.trace_inner_product(y
)
382 def natural_basis(self
):
384 Return a more-natural representation of this algebra's basis.
386 Every finite-dimensional Euclidean Jordan Algebra is a direct
387 sum of five simple algebras, four of which comprise Hermitian
388 matrices. This method returns the original "natural" basis
389 for our underlying vector space. (Typically, the natural basis
390 is used to construct the multiplication table in the first place.)
392 Note that this will always return a matrix. The standard basis
393 in `R^n` will be returned as `n`-by-`1` column matrices.
397 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
398 ....: RealSymmetricEJA)
402 sage: J = RealSymmetricEJA(2)
404 Finite family {0: e0, 1: e1, 2: e2}
405 sage: J.natural_basis()
413 sage: J = JordanSpinEJA(2)
415 Finite family {0: e0, 1: e1}
416 sage: J.natural_basis()
423 if self
._natural
_basis
is None:
424 return tuple( b
.to_vector().column() for b
in self
.basis() )
426 return self
._natural
_basis
432 Return the unit element of this algebra.
436 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
441 sage: J = RealCartesianProductEJA(5)
443 e0 + e1 + e2 + e3 + e4
447 The identity element acts like the identity::
449 sage: set_random_seed()
450 sage: J = random_eja()
451 sage: x = J.random_element()
452 sage: J.one()*x == x and x*J.one() == x
455 The matrix of the unit element's operator is the identity::
457 sage: set_random_seed()
458 sage: J = random_eja()
459 sage: actual = J.one().operator().matrix()
460 sage: expected = matrix.identity(J.base_ring(), J.dimension())
461 sage: actual == expected
465 # We can brute-force compute the matrices of the operators
466 # that correspond to the basis elements of this algebra.
467 # If some linear combination of those basis elements is the
468 # algebra identity, then the same linear combination of
469 # their matrices has to be the identity matrix.
471 # Of course, matrices aren't vectors in sage, so we have to
472 # appeal to the "long vectors" isometry.
473 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
475 # Now we use basis linear algebra to find the coefficients,
476 # of the matrices-as-vectors-linear-combination, which should
477 # work for the original algebra basis too.
478 A
= matrix
.column(self
.base_ring(), oper_vecs
)
480 # We used the isometry on the left-hand side already, but we
481 # still need to do it for the right-hand side. Recall that we
482 # wanted something that summed to the identity matrix.
483 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
485 # Now if there's an identity element in the algebra, this should work.
486 coeffs
= A
.solve_right(b
)
487 return self
.linear_combination(zip(self
.gens(), coeffs
))
492 Return the rank of this EJA.
496 The author knows of no algorithm to compute the rank of an EJA
497 where only the multiplication table is known. In lieu of one, we
498 require the rank to be specified when the algebra is created,
499 and simply pass along that number here.
503 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
504 ....: RealSymmetricEJA,
505 ....: ComplexHermitianEJA,
506 ....: QuaternionHermitianEJA,
511 The rank of the Jordan spin algebra is always two::
513 sage: JordanSpinEJA(2).rank()
515 sage: JordanSpinEJA(3).rank()
517 sage: JordanSpinEJA(4).rank()
520 The rank of the `n`-by-`n` Hermitian real, complex, or
521 quaternion matrices is `n`::
523 sage: RealSymmetricEJA(2).rank()
525 sage: ComplexHermitianEJA(2).rank()
527 sage: QuaternionHermitianEJA(2).rank()
529 sage: RealSymmetricEJA(5).rank()
531 sage: ComplexHermitianEJA(5).rank()
533 sage: QuaternionHermitianEJA(5).rank()
538 Ensure that every EJA that we know how to construct has a
539 positive integer rank::
541 sage: set_random_seed()
542 sage: r = random_eja().rank()
543 sage: r in ZZ and r > 0
550 def vector_space(self
):
552 Return the vector space that underlies this algebra.
556 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
560 sage: J = RealSymmetricEJA(2)
561 sage: J.vector_space()
562 Vector space of dimension 3 over Rational Field
565 return self
.zero().to_vector().parent().ambient_vector_space()
568 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
571 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
573 Return the Euclidean Jordan Algebra corresponding to the set
574 `R^n` under the Hadamard product.
576 Note: this is nothing more than the Cartesian product of ``n``
577 copies of the spin algebra. Once Cartesian product algebras
578 are implemented, this can go.
582 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
586 This multiplication table can be verified by hand::
588 sage: J = RealCartesianProductEJA(3)
589 sage: e0,e1,e2 = J.gens()
604 def __init__(self
, n
, field
=QQ
):
605 V
= VectorSpace(field
, n
)
606 mult_table
= [ [ V
.basis()[i
]*(i
== j
) for i
in range(n
) ]
609 fdeja
= super(RealCartesianProductEJA
, self
)
610 return fdeja
.__init
__(field
, mult_table
, rank
=n
)
612 def inner_product(self
, x
, y
):
613 return _usual_ip(x
,y
)
618 Return a "random" finite-dimensional Euclidean Jordan Algebra.
622 For now, we choose a random natural number ``n`` (greater than zero)
623 and then give you back one of the following:
625 * The cartesian product of the rational numbers ``n`` times; this is
626 ``QQ^n`` with the Hadamard product.
628 * The Jordan spin algebra on ``QQ^n``.
630 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
633 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
634 in the space of ``2n``-by-``2n`` real symmetric matrices.
636 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
637 in the space of ``4n``-by-``4n`` real symmetric matrices.
639 Later this might be extended to return Cartesian products of the
644 sage: from mjo.eja.eja_algebra import random_eja
649 Euclidean Jordan algebra of degree...
653 # The max_n component lets us choose different upper bounds on the
654 # value "n" that gets passed to the constructor. This is needed
655 # because e.g. R^{10} is reasonable to test, while the Hermitian
656 # 10-by-10 quaternion matrices are not.
657 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
659 (RealSymmetricEJA
, 5),
660 (ComplexHermitianEJA
, 4),
661 (QuaternionHermitianEJA
, 3)])
662 n
= ZZ
.random_element(1, max_n
)
663 return constructor(n
, field
=QQ
)
667 def _real_symmetric_basis(n
, field
=QQ
):
669 Return a basis for the space of real symmetric n-by-n matrices.
671 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
675 for j
in xrange(i
+1):
676 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
680 # Beware, orthogonal but not normalized!
681 Sij
= Eij
+ Eij
.transpose()
686 def _complex_hermitian_basis(n
, field
=QQ
):
688 Returns a basis for the space of complex Hermitian n-by-n matrices.
692 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
696 sage: set_random_seed()
697 sage: n = ZZ.random_element(1,5)
698 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
702 F
= QuadraticField(-1, 'I')
705 # This is like the symmetric case, but we need to be careful:
707 # * We want conjugate-symmetry, not just symmetry.
708 # * The diagonal will (as a result) be real.
712 for j
in xrange(i
+1):
713 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
715 Sij
= _embed_complex_matrix(Eij
)
718 # Beware, orthogonal but not normalized! The second one
719 # has a minus because it's conjugated.
720 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
722 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
727 def _quaternion_hermitian_basis(n
, field
=QQ
):
729 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
733 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
737 sage: set_random_seed()
738 sage: n = ZZ.random_element(1,5)
739 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
743 Q
= QuaternionAlgebra(QQ
,-1,-1)
746 # This is like the symmetric case, but we need to be careful:
748 # * We want conjugate-symmetry, not just symmetry.
749 # * The diagonal will (as a result) be real.
753 for j
in xrange(i
+1):
754 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
756 Sij
= _embed_quaternion_matrix(Eij
)
759 # Beware, orthogonal but not normalized! The second,
760 # third, and fourth ones have a minus because they're
762 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
764 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
766 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
768 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
774 def _multiplication_table_from_matrix_basis(basis
):
776 At least three of the five simple Euclidean Jordan algebras have the
777 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
778 multiplication on the right is matrix multiplication. Given a basis
779 for the underlying matrix space, this function returns a
780 multiplication table (obtained by looping through the basis
781 elements) for an algebra of those matrices.
783 # In S^2, for example, we nominally have four coordinates even
784 # though the space is of dimension three only. The vector space V
785 # is supposed to hold the entire long vector, and the subspace W
786 # of V will be spanned by the vectors that arise from symmetric
787 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
788 field
= basis
[0].base_ring()
789 dimension
= basis
[0].nrows()
791 V
= VectorSpace(field
, dimension
**2)
792 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
794 mult_table
= [[W
.zero() for i
in range(n
)] for j
in range(n
)]
797 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
798 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
803 def _embed_complex_matrix(M
):
805 Embed the n-by-n complex matrix ``M`` into the space of real
806 matrices of size 2n-by-2n via the map the sends each entry `z = a +
807 bi` to the block matrix ``[[a,b],[-b,a]]``.
811 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
815 sage: F = QuadraticField(-1,'i')
816 sage: x1 = F(4 - 2*i)
817 sage: x2 = F(1 + 2*i)
820 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
821 sage: _embed_complex_matrix(M)
830 Embedding is a homomorphism (isomorphism, in fact)::
832 sage: set_random_seed()
833 sage: n = ZZ.random_element(5)
834 sage: F = QuadraticField(-1, 'i')
835 sage: X = random_matrix(F, n)
836 sage: Y = random_matrix(F, n)
837 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
838 sage: expected = _embed_complex_matrix(X*Y)
839 sage: actual == expected
845 raise ValueError("the matrix 'M' must be square")
846 field
= M
.base_ring()
851 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
853 # We can drop the imaginaries here.
854 return matrix
.block(field
.base_ring(), n
, blocks
)
857 def _unembed_complex_matrix(M
):
859 The inverse of _embed_complex_matrix().
863 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
864 ....: _unembed_complex_matrix)
868 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
869 ....: [-2, 1, -4, 3],
870 ....: [ 9, 10, 11, 12],
871 ....: [-10, 9, -12, 11] ])
872 sage: _unembed_complex_matrix(A)
874 [ 10*i + 9 12*i + 11]
878 Unembedding is the inverse of embedding::
880 sage: set_random_seed()
881 sage: F = QuadraticField(-1, 'i')
882 sage: M = random_matrix(F, 3)
883 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
889 raise ValueError("the matrix 'M' must be square")
890 if not n
.mod(2).is_zero():
891 raise ValueError("the matrix 'M' must be a complex embedding")
893 F
= QuadraticField(-1, 'i')
896 # Go top-left to bottom-right (reading order), converting every
897 # 2-by-2 block we see to a single complex element.
899 for k
in xrange(n
/2):
900 for j
in xrange(n
/2):
901 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
902 if submat
[0,0] != submat
[1,1]:
903 raise ValueError('bad on-diagonal submatrix')
904 if submat
[0,1] != -submat
[1,0]:
905 raise ValueError('bad off-diagonal submatrix')
906 z
= submat
[0,0] + submat
[0,1]*i
909 return matrix(F
, n
/2, elements
)
912 def _embed_quaternion_matrix(M
):
914 Embed the n-by-n quaternion matrix ``M`` into the space of real
915 matrices of size 4n-by-4n by first sending each quaternion entry
916 `z = a + bi + cj + dk` to the block-complex matrix
917 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
922 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
926 sage: Q = QuaternionAlgebra(QQ,-1,-1)
927 sage: i,j,k = Q.gens()
928 sage: x = 1 + 2*i + 3*j + 4*k
929 sage: M = matrix(Q, 1, [[x]])
930 sage: _embed_quaternion_matrix(M)
936 Embedding is a homomorphism (isomorphism, in fact)::
938 sage: set_random_seed()
939 sage: n = ZZ.random_element(5)
940 sage: Q = QuaternionAlgebra(QQ,-1,-1)
941 sage: X = random_matrix(Q, n)
942 sage: Y = random_matrix(Q, n)
943 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
944 sage: expected = _embed_quaternion_matrix(X*Y)
945 sage: actual == expected
949 quaternions
= M
.base_ring()
952 raise ValueError("the matrix 'M' must be square")
954 F
= QuadraticField(-1, 'i')
959 t
= z
.coefficient_tuple()
964 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
965 [-c
+ d
*i
, a
- b
*i
]])
966 blocks
.append(_embed_complex_matrix(cplx_matrix
))
968 # We should have real entries by now, so use the realest field
969 # we've got for the return value.
970 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
973 def _unembed_quaternion_matrix(M
):
975 The inverse of _embed_quaternion_matrix().
979 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
980 ....: _unembed_quaternion_matrix)
984 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
985 ....: [-2, 1, -4, 3],
986 ....: [-3, 4, 1, -2],
987 ....: [-4, -3, 2, 1]])
988 sage: _unembed_quaternion_matrix(M)
989 [1 + 2*i + 3*j + 4*k]
993 Unembedding is the inverse of embedding::
995 sage: set_random_seed()
996 sage: Q = QuaternionAlgebra(QQ, -1, -1)
997 sage: M = random_matrix(Q, 3)
998 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1004 raise ValueError("the matrix 'M' must be square")
1005 if not n
.mod(4).is_zero():
1006 raise ValueError("the matrix 'M' must be a complex embedding")
1008 Q
= QuaternionAlgebra(QQ
,-1,-1)
1011 # Go top-left to bottom-right (reading order), converting every
1012 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1015 for l
in xrange(n
/4):
1016 for m
in xrange(n
/4):
1017 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1018 if submat
[0,0] != submat
[1,1].conjugate():
1019 raise ValueError('bad on-diagonal submatrix')
1020 if submat
[0,1] != -submat
[1,0].conjugate():
1021 raise ValueError('bad off-diagonal submatrix')
1022 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1023 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1026 return matrix(Q
, n
/4, elements
)
1029 # The usual inner product on R^n.
1031 return x
.to_vector().inner_product(y
.to_vector())
1033 # The inner product used for the real symmetric simple EJA.
1034 # We keep it as a separate function because e.g. the complex
1035 # algebra uses the same inner product, except divided by 2.
1036 def _matrix_ip(X
,Y
):
1037 X_mat
= X
.natural_representation()
1038 Y_mat
= Y
.natural_representation()
1039 return (X_mat
*Y_mat
).trace()
1042 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1044 The rank-n simple EJA consisting of real symmetric n-by-n
1045 matrices, the usual symmetric Jordan product, and the trace inner
1046 product. It has dimension `(n^2 + n)/2` over the reals.
1050 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1054 sage: J = RealSymmetricEJA(2)
1055 sage: e0, e1, e2 = J.gens()
1065 The dimension of this algebra is `(n^2 + n) / 2`::
1067 sage: set_random_seed()
1068 sage: n = ZZ.random_element(1,5)
1069 sage: J = RealSymmetricEJA(n)
1070 sage: J.dimension() == (n^2 + n)/2
1073 The Jordan multiplication is what we think it is::
1075 sage: set_random_seed()
1076 sage: n = ZZ.random_element(1,5)
1077 sage: J = RealSymmetricEJA(n)
1078 sage: x = J.random_element()
1079 sage: y = J.random_element()
1080 sage: actual = (x*y).natural_representation()
1081 sage: X = x.natural_representation()
1082 sage: Y = y.natural_representation()
1083 sage: expected = (X*Y + Y*X)/2
1084 sage: actual == expected
1086 sage: J(expected) == x*y
1090 def __init__(self
, n
, field
=QQ
):
1091 S
= _real_symmetric_basis(n
, field
=field
)
1092 Qs
= _multiplication_table_from_matrix_basis(S
)
1094 fdeja
= super(RealSymmetricEJA
, self
)
1095 return fdeja
.__init
__(field
,
1100 def inner_product(self
, x
, y
):
1101 return _matrix_ip(x
,y
)
1104 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1106 The rank-n simple EJA consisting of complex Hermitian n-by-n
1107 matrices over the real numbers, the usual symmetric Jordan product,
1108 and the real-part-of-trace inner product. It has dimension `n^2` over
1113 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1117 The dimension of this algebra is `n^2`::
1119 sage: set_random_seed()
1120 sage: n = ZZ.random_element(1,5)
1121 sage: J = ComplexHermitianEJA(n)
1122 sage: J.dimension() == n^2
1125 The Jordan multiplication is what we think it is::
1127 sage: set_random_seed()
1128 sage: n = ZZ.random_element(1,5)
1129 sage: J = ComplexHermitianEJA(n)
1130 sage: x = J.random_element()
1131 sage: y = J.random_element()
1132 sage: actual = (x*y).natural_representation()
1133 sage: X = x.natural_representation()
1134 sage: Y = y.natural_representation()
1135 sage: expected = (X*Y + Y*X)/2
1136 sage: actual == expected
1138 sage: J(expected) == x*y
1142 def __init__(self
, n
, field
=QQ
):
1143 S
= _complex_hermitian_basis(n
)
1144 Qs
= _multiplication_table_from_matrix_basis(S
)
1146 fdeja
= super(ComplexHermitianEJA
, self
)
1147 return fdeja
.__init
__(field
,
1153 def inner_product(self
, x
, y
):
1154 # Since a+bi on the diagonal is represented as
1159 # we'll double-count the "a" entries if we take the trace of
1161 return _matrix_ip(x
,y
)/2
1164 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1166 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1167 matrices, the usual symmetric Jordan product, and the
1168 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1173 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1177 The dimension of this algebra is `n^2`::
1179 sage: set_random_seed()
1180 sage: n = ZZ.random_element(1,5)
1181 sage: J = QuaternionHermitianEJA(n)
1182 sage: J.dimension() == 2*(n^2) - n
1185 The Jordan multiplication is what we think it is::
1187 sage: set_random_seed()
1188 sage: n = ZZ.random_element(1,5)
1189 sage: J = QuaternionHermitianEJA(n)
1190 sage: x = J.random_element()
1191 sage: y = J.random_element()
1192 sage: actual = (x*y).natural_representation()
1193 sage: X = x.natural_representation()
1194 sage: Y = y.natural_representation()
1195 sage: expected = (X*Y + Y*X)/2
1196 sage: actual == expected
1198 sage: J(expected) == x*y
1202 def __init__(self
, n
, field
=QQ
):
1203 S
= _quaternion_hermitian_basis(n
)
1204 Qs
= _multiplication_table_from_matrix_basis(S
)
1206 fdeja
= super(QuaternionHermitianEJA
, self
)
1207 return fdeja
.__init
__(field
,
1212 def inner_product(self
, x
, y
):
1213 # Since a+bi+cj+dk on the diagonal is represented as
1215 # a + bi +cj + dk = [ a b c d]
1220 # we'll quadruple-count the "a" entries if we take the trace of
1222 return _matrix_ip(x
,y
)/4
1225 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1227 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1228 with the usual inner product and jordan product ``x*y =
1229 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1234 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1238 This multiplication table can be verified by hand::
1240 sage: J = JordanSpinEJA(4)
1241 sage: e0,e1,e2,e3 = J.gens()
1258 def __init__(self
, n
, field
=QQ
):
1259 V
= VectorSpace(field
, n
)
1260 mult_table
= [[V
.zero() for i
in range(n
)] for j
in range(n
)]
1270 z0
= x
.inner_product(y
)
1271 zbar
= y0
*xbar
+ x0
*ybar
1272 z
= V([z0
] + zbar
.list())
1273 mult_table
[i
][j
] = z
1275 # The rank of the spin algebra is two, unless we're in a
1276 # one-dimensional ambient space (because the rank is bounded by
1277 # the ambient dimension).
1278 fdeja
= super(JordanSpinEJA
, self
)
1279 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2))
1281 def inner_product(self
, x
, y
):
1282 return _usual_ip(x
,y
)