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
= matrix(
68 [ map(lambda x
: self
.from_vector(x
), ls
)
69 for ls
in mult_table
] )
70 self
._multiplication
_table
.set_immutable()
73 def _element_constructor_(self
, elt
):
75 Construct an element of this algebra from its natural
78 This gets called only after the parent element _call_ method
79 fails to find a coercion for the argument.
83 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
84 ....: RealCartesianProductEJA,
85 ....: RealSymmetricEJA)
89 The identity in `S^n` is converted to the identity in the EJA::
91 sage: J = RealSymmetricEJA(3)
92 sage: I = matrix.identity(QQ,3)
96 This skew-symmetric matrix can't be represented in the EJA::
98 sage: J = RealSymmetricEJA(3)
99 sage: A = matrix(QQ,3, lambda i,j: i-j)
101 Traceback (most recent call last):
103 ArithmeticError: vector is not in free module
107 Ensure that we can convert any element of the two non-matrix
108 simple algebras (whose natural representations are their usual
109 vector representations) back and forth faithfully::
111 sage: set_random_seed()
112 sage: J = RealCartesianProductEJA(5)
113 sage: x = J.random_element()
114 sage: J(x.to_vector().column()) == x
116 sage: J = JordanSpinEJA(5)
117 sage: x = J.random_element()
118 sage: J(x.to_vector().column()) == x
122 natural_basis
= self
.natural_basis()
123 if elt
not in natural_basis
[0].matrix_space():
124 raise ValueError("not a naturally-represented algebra element")
126 # Thanks for nothing! Matrix spaces aren't vector
127 # spaces in Sage, so we have to figure out its
128 # natural-basis coordinates ourselves.
129 V
= VectorSpace(elt
.base_ring(), elt
.nrows()*elt
.ncols())
130 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
131 coords
= W
.coordinate_vector(_mat2vec(elt
))
132 return self
.from_vector(coords
)
137 Return a string representation of ``self``.
141 sage: from mjo.eja.eja_algebra import JordanSpinEJA
145 Ensure that it says what we think it says::
147 sage: JordanSpinEJA(2, field=QQ)
148 Euclidean Jordan algebra of degree 2 over Rational Field
149 sage: JordanSpinEJA(3, field=RDF)
150 Euclidean Jordan algebra of degree 3 over Real Double Field
153 # TODO: change this to say "dimension" and fix all the tests.
154 fmt
= "Euclidean Jordan algebra of degree {} over {}"
155 return fmt
.format(self
.dimension(), self
.base_ring())
157 def product_on_basis(self
, i
, j
):
158 return self
._multiplication
_table
[i
,j
]
160 def _a_regular_element(self
):
162 Guess a regular element. Needed to compute the basis for our
163 characteristic polynomial coefficients.
167 sage: from mjo.eja.eja_algebra import random_eja
171 Ensure that this hacky method succeeds for every algebra that we
172 know how to construct::
174 sage: set_random_seed()
175 sage: J = random_eja()
176 sage: J._a_regular_element().is_regular()
181 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
182 if not z
.is_regular():
183 raise ValueError("don't know a regular element")
188 def _charpoly_basis_space(self
):
190 Return the vector space spanned by the basis used in our
191 characteristic polynomial coefficients. This is used not only to
192 compute those coefficients, but also any time we need to
193 evaluate the coefficients (like when we compute the trace or
196 z
= self
._a
_regular
_element
()
197 V
= self
.vector_space()
198 V1
= V
.span_of_basis( (z
**k
).to_vector() for k
in range(self
.rank()) )
199 b
= (V1
.basis() + V1
.complement().basis())
200 return V
.span_of_basis(b
)
204 def _charpoly_coeff(self
, i
):
206 Return the coefficient polynomial "a_{i}" of this algebra's
207 general characteristic polynomial.
209 Having this be a separate cached method lets us compute and
210 store the trace/determinant (a_{r-1} and a_{0} respectively)
211 separate from the entire characteristic polynomial.
213 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
214 R
= A_of_x
.base_ring()
216 # Guaranteed by theory
219 # Danger: the in-place modification is done for performance
220 # reasons (reconstructing a matrix with huge polynomial
221 # entries is slow), but I don't know how cached_method works,
222 # so it's highly possible that we're modifying some global
223 # list variable by reference, here. In other words, you
224 # probably shouldn't call this method twice on the same
225 # algebra, at the same time, in two threads
226 Ai_orig
= A_of_x
.column(i
)
227 A_of_x
.set_column(i
,xr
)
228 numerator
= A_of_x
.det()
229 A_of_x
.set_column(i
,Ai_orig
)
231 # We're relying on the theory here to ensure that each a_i is
232 # indeed back in R, and the added negative signs are to make
233 # the whole charpoly expression sum to zero.
234 return R(-numerator
/detA
)
238 def _charpoly_matrix_system(self
):
240 Compute the matrix whose entries A_ij are polynomials in
241 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
242 corresponding to `x^r` and the determinent of the matrix A =
243 [A_ij]. In other words, all of the fixed (cachable) data needed
244 to compute the coefficients of the characteristic polynomial.
249 # Construct a new algebra over a multivariate polynomial ring...
250 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
251 R
= PolynomialRing(self
.base_ring(), names
)
252 # Hack around the fact that our multiplication table is in terms of
253 # algebra elements but the constructor wants it in terms of vectors.
254 vmt
= [ tuple([ self
._multiplication
_table
[i
,j
].to_vector()
255 for j
in range(self
._multiplication
_table
.nrows()) ])
256 for i
in range(self
._multiplication
_table
.ncols()) ]
257 J
= FiniteDimensionalEuclideanJordanAlgebra(R
, tuple(vmt
), r
)
259 idmat
= matrix
.identity(J
.base_ring(), n
)
261 W
= self
._charpoly
_basis
_space
()
262 W
= W
.change_ring(R
.fraction_field())
264 # Starting with the standard coordinates x = (X1,X2,...,Xn)
265 # and then converting the entries to W-coordinates allows us
266 # to pass in the standard coordinates to the charpoly and get
267 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
270 # W.coordinates(x^2) eval'd at (standard z-coords)
274 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
276 # We want the middle equivalent thing in our matrix, but use
277 # the first equivalent thing instead so that we can pass in
278 # standard coordinates.
279 x
= J
.from_vector(W(R
.gens()))
281 # Handle the zeroth power separately, because computing
282 # the unit element in J is mathematically suspect.
283 x0
= W
.coordinate_vector(self
.one().to_vector())
285 l1
+= [ W
.coordinate_vector((x
**k
).to_vector()).column()
286 for k
in range(1,r
) ]
287 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
288 A_of_x
= matrix
.block(R
, 1, n
, (l1
+ l2
))
289 xr
= W
.coordinate_vector((x
**r
).to_vector())
290 return (A_of_x
, x
, xr
, A_of_x
.det())
294 def characteristic_polynomial(self
):
296 Return a characteristic polynomial that works for all elements
299 The resulting polynomial has `n+1` variables, where `n` is the
300 dimension of this algebra. The first `n` variables correspond to
301 the coordinates of an algebra element: when evaluated at the
302 coordinates of an algebra element with respect to a certain
303 basis, the result is a univariate polynomial (in the one
304 remaining variable ``t``), namely the characteristic polynomial
309 sage: from mjo.eja.eja_algebra import JordanSpinEJA
313 The characteristic polynomial in the spin algebra is given in
314 Alizadeh, Example 11.11::
316 sage: J = JordanSpinEJA(3)
317 sage: p = J.characteristic_polynomial(); p
318 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
319 sage: xvec = J.one().to_vector()
327 # The list of coefficient polynomials a_1, a_2, ..., a_n.
328 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
330 # We go to a bit of trouble here to reorder the
331 # indeterminates, so that it's easier to evaluate the
332 # characteristic polynomial at x's coordinates and get back
333 # something in terms of t, which is what we want.
335 S
= PolynomialRing(self
.base_ring(),'t')
337 S
= PolynomialRing(S
, R
.variable_names())
340 # Note: all entries past the rth should be zero. The
341 # coefficient of the highest power (x^r) is 1, but it doesn't
342 # appear in the solution vector which contains coefficients
343 # for the other powers (to make them sum to x^r).
345 a
[r
] = 1 # corresponds to x^r
347 # When the rank is equal to the dimension, trying to
348 # assign a[r] goes out-of-bounds.
349 a
.append(1) # corresponds to x^r
351 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
354 def inner_product(self
, x
, y
):
356 The inner product associated with this Euclidean Jordan algebra.
358 Defaults to the trace inner product, but can be overridden by
359 subclasses if they are sure that the necessary properties are
364 sage: from mjo.eja.eja_algebra import random_eja
368 The inner product must satisfy its axiom for this algebra to truly
369 be a Euclidean Jordan Algebra::
371 sage: set_random_seed()
372 sage: J = random_eja()
373 sage: x = J.random_element()
374 sage: y = J.random_element()
375 sage: z = J.random_element()
376 sage: (x*y).inner_product(z) == y.inner_product(x*z)
380 if (not x
in self
) or (not y
in self
):
381 raise TypeError("arguments must live in this algebra")
382 return x
.trace_inner_product(y
)
385 def multiplication_table(self
):
387 Return a readable matrix representation of this algebra's
388 multiplication table. The (i,j)th entry in the matrix contains
389 the product of the ith basis element with the jth.
391 This is not extraordinarily useful, but it overrides a superclass
392 method that would otherwise just crash and complain about the
393 algebra being infinite.
397 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
398 ....: RealCartesianProductEJA)
402 sage: J = RealCartesianProductEJA(3)
403 sage: J.multiplication_table()
410 sage: J = JordanSpinEJA(3)
411 sage: J.multiplication_table()
417 return self
._multiplication
_table
420 def natural_basis(self
):
422 Return a more-natural representation of this algebra's basis.
424 Every finite-dimensional Euclidean Jordan Algebra is a direct
425 sum of five simple algebras, four of which comprise Hermitian
426 matrices. This method returns the original "natural" basis
427 for our underlying vector space. (Typically, the natural basis
428 is used to construct the multiplication table in the first place.)
430 Note that this will always return a matrix. The standard basis
431 in `R^n` will be returned as `n`-by-`1` column matrices.
435 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
436 ....: RealSymmetricEJA)
440 sage: J = RealSymmetricEJA(2)
442 Finite family {0: e0, 1: e1, 2: e2}
443 sage: J.natural_basis()
451 sage: J = JordanSpinEJA(2)
453 Finite family {0: e0, 1: e1}
454 sage: J.natural_basis()
461 if self
._natural
_basis
is None:
462 return tuple( b
.to_vector().column() for b
in self
.basis() )
464 return self
._natural
_basis
470 Return the unit element of this algebra.
474 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
479 sage: J = RealCartesianProductEJA(5)
481 e0 + e1 + e2 + e3 + e4
485 The identity element acts like the identity::
487 sage: set_random_seed()
488 sage: J = random_eja()
489 sage: x = J.random_element()
490 sage: J.one()*x == x and x*J.one() == x
493 The matrix of the unit element's operator is the identity::
495 sage: set_random_seed()
496 sage: J = random_eja()
497 sage: actual = J.one().operator().matrix()
498 sage: expected = matrix.identity(J.base_ring(), J.dimension())
499 sage: actual == expected
503 # We can brute-force compute the matrices of the operators
504 # that correspond to the basis elements of this algebra.
505 # If some linear combination of those basis elements is the
506 # algebra identity, then the same linear combination of
507 # their matrices has to be the identity matrix.
509 # Of course, matrices aren't vectors in sage, so we have to
510 # appeal to the "long vectors" isometry.
511 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
513 # Now we use basis linear algebra to find the coefficients,
514 # of the matrices-as-vectors-linear-combination, which should
515 # work for the original algebra basis too.
516 A
= matrix
.column(self
.base_ring(), oper_vecs
)
518 # We used the isometry on the left-hand side already, but we
519 # still need to do it for the right-hand side. Recall that we
520 # wanted something that summed to the identity matrix.
521 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
523 # Now if there's an identity element in the algebra, this should work.
524 coeffs
= A
.solve_right(b
)
525 return self
.linear_combination(zip(self
.gens(), coeffs
))
530 Return the rank of this EJA.
534 The author knows of no algorithm to compute the rank of an EJA
535 where only the multiplication table is known. In lieu of one, we
536 require the rank to be specified when the algebra is created,
537 and simply pass along that number here.
541 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
542 ....: RealSymmetricEJA,
543 ....: ComplexHermitianEJA,
544 ....: QuaternionHermitianEJA,
549 The rank of the Jordan spin algebra is always two::
551 sage: JordanSpinEJA(2).rank()
553 sage: JordanSpinEJA(3).rank()
555 sage: JordanSpinEJA(4).rank()
558 The rank of the `n`-by-`n` Hermitian real, complex, or
559 quaternion matrices is `n`::
561 sage: RealSymmetricEJA(2).rank()
563 sage: ComplexHermitianEJA(2).rank()
565 sage: QuaternionHermitianEJA(2).rank()
567 sage: RealSymmetricEJA(5).rank()
569 sage: ComplexHermitianEJA(5).rank()
571 sage: QuaternionHermitianEJA(5).rank()
576 Ensure that every EJA that we know how to construct has a
577 positive integer rank::
579 sage: set_random_seed()
580 sage: r = random_eja().rank()
581 sage: r in ZZ and r > 0
588 def vector_space(self
):
590 Return the vector space that underlies this algebra.
594 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
598 sage: J = RealSymmetricEJA(2)
599 sage: J.vector_space()
600 Vector space of dimension 3 over Rational Field
603 return self
.zero().to_vector().parent().ambient_vector_space()
606 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
609 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
611 Return the Euclidean Jordan Algebra corresponding to the set
612 `R^n` under the Hadamard product.
614 Note: this is nothing more than the Cartesian product of ``n``
615 copies of the spin algebra. Once Cartesian product algebras
616 are implemented, this can go.
620 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
624 This multiplication table can be verified by hand::
626 sage: J = RealCartesianProductEJA(3)
627 sage: e0,e1,e2 = J.gens()
642 def __init__(self
, n
, field
=QQ
):
643 V
= VectorSpace(field
, n
)
644 mult_table
= [ [ V
.basis()[i
]*(i
== j
) for j
in range(n
) ]
647 fdeja
= super(RealCartesianProductEJA
, self
)
648 return fdeja
.__init
__(field
, mult_table
, rank
=n
)
650 def inner_product(self
, x
, y
):
651 return _usual_ip(x
,y
)
656 Return a "random" finite-dimensional Euclidean Jordan Algebra.
660 For now, we choose a random natural number ``n`` (greater than zero)
661 and then give you back one of the following:
663 * The cartesian product of the rational numbers ``n`` times; this is
664 ``QQ^n`` with the Hadamard product.
666 * The Jordan spin algebra on ``QQ^n``.
668 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
671 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
672 in the space of ``2n``-by-``2n`` real symmetric matrices.
674 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
675 in the space of ``4n``-by-``4n`` real symmetric matrices.
677 Later this might be extended to return Cartesian products of the
682 sage: from mjo.eja.eja_algebra import random_eja
687 Euclidean Jordan algebra of degree...
691 # The max_n component lets us choose different upper bounds on the
692 # value "n" that gets passed to the constructor. This is needed
693 # because e.g. R^{10} is reasonable to test, while the Hermitian
694 # 10-by-10 quaternion matrices are not.
695 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
697 (RealSymmetricEJA
, 5),
698 (ComplexHermitianEJA
, 4),
699 (QuaternionHermitianEJA
, 3)])
700 n
= ZZ
.random_element(1, max_n
)
701 return constructor(n
, field
=QQ
)
705 def _real_symmetric_basis(n
, field
=QQ
):
707 Return a basis for the space of real symmetric n-by-n matrices.
709 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
713 for j
in xrange(i
+1):
714 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
718 # Beware, orthogonal but not normalized!
719 Sij
= Eij
+ Eij
.transpose()
724 def _complex_hermitian_basis(n
, field
=QQ
):
726 Returns a basis for the space of complex Hermitian n-by-n matrices.
730 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
734 sage: set_random_seed()
735 sage: n = ZZ.random_element(1,5)
736 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
740 F
= QuadraticField(-1, 'I')
743 # This is like the symmetric case, but we need to be careful:
745 # * We want conjugate-symmetry, not just symmetry.
746 # * The diagonal will (as a result) be real.
750 for j
in xrange(i
+1):
751 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
753 Sij
= _embed_complex_matrix(Eij
)
756 # Beware, orthogonal but not normalized! The second one
757 # has a minus because it's conjugated.
758 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
760 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
765 def _quaternion_hermitian_basis(n
, field
=QQ
):
767 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
771 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
775 sage: set_random_seed()
776 sage: n = ZZ.random_element(1,5)
777 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
781 Q
= QuaternionAlgebra(QQ
,-1,-1)
784 # This is like the symmetric case, but we need to be careful:
786 # * We want conjugate-symmetry, not just symmetry.
787 # * The diagonal will (as a result) be real.
791 for j
in xrange(i
+1):
792 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
794 Sij
= _embed_quaternion_matrix(Eij
)
797 # Beware, orthogonal but not normalized! The second,
798 # third, and fourth ones have a minus because they're
800 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
802 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
804 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
806 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
812 def _multiplication_table_from_matrix_basis(basis
):
814 At least three of the five simple Euclidean Jordan algebras have the
815 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
816 multiplication on the right is matrix multiplication. Given a basis
817 for the underlying matrix space, this function returns a
818 multiplication table (obtained by looping through the basis
819 elements) for an algebra of those matrices.
821 # In S^2, for example, we nominally have four coordinates even
822 # though the space is of dimension three only. The vector space V
823 # is supposed to hold the entire long vector, and the subspace W
824 # of V will be spanned by the vectors that arise from symmetric
825 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
826 field
= basis
[0].base_ring()
827 dimension
= basis
[0].nrows()
829 V
= VectorSpace(field
, dimension
**2)
830 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
832 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
835 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
836 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
841 def _embed_complex_matrix(M
):
843 Embed the n-by-n complex matrix ``M`` into the space of real
844 matrices of size 2n-by-2n via the map the sends each entry `z = a +
845 bi` to the block matrix ``[[a,b],[-b,a]]``.
849 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
853 sage: F = QuadraticField(-1,'i')
854 sage: x1 = F(4 - 2*i)
855 sage: x2 = F(1 + 2*i)
858 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
859 sage: _embed_complex_matrix(M)
868 Embedding is a homomorphism (isomorphism, in fact)::
870 sage: set_random_seed()
871 sage: n = ZZ.random_element(5)
872 sage: F = QuadraticField(-1, 'i')
873 sage: X = random_matrix(F, n)
874 sage: Y = random_matrix(F, n)
875 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
876 sage: expected = _embed_complex_matrix(X*Y)
877 sage: actual == expected
883 raise ValueError("the matrix 'M' must be square")
884 field
= M
.base_ring()
889 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
891 # We can drop the imaginaries here.
892 return matrix
.block(field
.base_ring(), n
, blocks
)
895 def _unembed_complex_matrix(M
):
897 The inverse of _embed_complex_matrix().
901 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
902 ....: _unembed_complex_matrix)
906 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
907 ....: [-2, 1, -4, 3],
908 ....: [ 9, 10, 11, 12],
909 ....: [-10, 9, -12, 11] ])
910 sage: _unembed_complex_matrix(A)
912 [ 10*i + 9 12*i + 11]
916 Unembedding is the inverse of embedding::
918 sage: set_random_seed()
919 sage: F = QuadraticField(-1, 'i')
920 sage: M = random_matrix(F, 3)
921 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
927 raise ValueError("the matrix 'M' must be square")
928 if not n
.mod(2).is_zero():
929 raise ValueError("the matrix 'M' must be a complex embedding")
931 F
= QuadraticField(-1, 'i')
934 # Go top-left to bottom-right (reading order), converting every
935 # 2-by-2 block we see to a single complex element.
937 for k
in xrange(n
/2):
938 for j
in xrange(n
/2):
939 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
940 if submat
[0,0] != submat
[1,1]:
941 raise ValueError('bad on-diagonal submatrix')
942 if submat
[0,1] != -submat
[1,0]:
943 raise ValueError('bad off-diagonal submatrix')
944 z
= submat
[0,0] + submat
[0,1]*i
947 return matrix(F
, n
/2, elements
)
950 def _embed_quaternion_matrix(M
):
952 Embed the n-by-n quaternion matrix ``M`` into the space of real
953 matrices of size 4n-by-4n by first sending each quaternion entry
954 `z = a + bi + cj + dk` to the block-complex matrix
955 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
960 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
964 sage: Q = QuaternionAlgebra(QQ,-1,-1)
965 sage: i,j,k = Q.gens()
966 sage: x = 1 + 2*i + 3*j + 4*k
967 sage: M = matrix(Q, 1, [[x]])
968 sage: _embed_quaternion_matrix(M)
974 Embedding is a homomorphism (isomorphism, in fact)::
976 sage: set_random_seed()
977 sage: n = ZZ.random_element(5)
978 sage: Q = QuaternionAlgebra(QQ,-1,-1)
979 sage: X = random_matrix(Q, n)
980 sage: Y = random_matrix(Q, n)
981 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
982 sage: expected = _embed_quaternion_matrix(X*Y)
983 sage: actual == expected
987 quaternions
= M
.base_ring()
990 raise ValueError("the matrix 'M' must be square")
992 F
= QuadraticField(-1, 'i')
997 t
= z
.coefficient_tuple()
1002 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1003 [-c
+ d
*i
, a
- b
*i
]])
1004 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1006 # We should have real entries by now, so use the realest field
1007 # we've got for the return value.
1008 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1011 def _unembed_quaternion_matrix(M
):
1013 The inverse of _embed_quaternion_matrix().
1017 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1018 ....: _unembed_quaternion_matrix)
1022 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1023 ....: [-2, 1, -4, 3],
1024 ....: [-3, 4, 1, -2],
1025 ....: [-4, -3, 2, 1]])
1026 sage: _unembed_quaternion_matrix(M)
1027 [1 + 2*i + 3*j + 4*k]
1031 Unembedding is the inverse of embedding::
1033 sage: set_random_seed()
1034 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1035 sage: M = random_matrix(Q, 3)
1036 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1042 raise ValueError("the matrix 'M' must be square")
1043 if not n
.mod(4).is_zero():
1044 raise ValueError("the matrix 'M' must be a complex embedding")
1046 Q
= QuaternionAlgebra(QQ
,-1,-1)
1049 # Go top-left to bottom-right (reading order), converting every
1050 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1053 for l
in xrange(n
/4):
1054 for m
in xrange(n
/4):
1055 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1056 if submat
[0,0] != submat
[1,1].conjugate():
1057 raise ValueError('bad on-diagonal submatrix')
1058 if submat
[0,1] != -submat
[1,0].conjugate():
1059 raise ValueError('bad off-diagonal submatrix')
1060 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1061 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1064 return matrix(Q
, n
/4, elements
)
1067 # The usual inner product on R^n.
1069 return x
.to_vector().inner_product(y
.to_vector())
1071 # The inner product used for the real symmetric simple EJA.
1072 # We keep it as a separate function because e.g. the complex
1073 # algebra uses the same inner product, except divided by 2.
1074 def _matrix_ip(X
,Y
):
1075 X_mat
= X
.natural_representation()
1076 Y_mat
= Y
.natural_representation()
1077 return (X_mat
*Y_mat
).trace()
1080 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1082 The rank-n simple EJA consisting of real symmetric n-by-n
1083 matrices, the usual symmetric Jordan product, and the trace inner
1084 product. It has dimension `(n^2 + n)/2` over the reals.
1088 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1092 sage: J = RealSymmetricEJA(2)
1093 sage: e0, e1, e2 = J.gens()
1103 The dimension of this algebra is `(n^2 + n) / 2`::
1105 sage: set_random_seed()
1106 sage: n = ZZ.random_element(1,5)
1107 sage: J = RealSymmetricEJA(n)
1108 sage: J.dimension() == (n^2 + n)/2
1111 The Jordan multiplication is what we think it is::
1113 sage: set_random_seed()
1114 sage: n = ZZ.random_element(1,5)
1115 sage: J = RealSymmetricEJA(n)
1116 sage: x = J.random_element()
1117 sage: y = J.random_element()
1118 sage: actual = (x*y).natural_representation()
1119 sage: X = x.natural_representation()
1120 sage: Y = y.natural_representation()
1121 sage: expected = (X*Y + Y*X)/2
1122 sage: actual == expected
1124 sage: J(expected) == x*y
1128 def __init__(self
, n
, field
=QQ
):
1129 S
= _real_symmetric_basis(n
, field
=field
)
1130 Qs
= _multiplication_table_from_matrix_basis(S
)
1132 fdeja
= super(RealSymmetricEJA
, self
)
1133 return fdeja
.__init
__(field
,
1138 def inner_product(self
, x
, y
):
1139 return _matrix_ip(x
,y
)
1142 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1144 The rank-n simple EJA consisting of complex Hermitian n-by-n
1145 matrices over the real numbers, the usual symmetric Jordan product,
1146 and the real-part-of-trace inner product. It has dimension `n^2` over
1151 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1155 The dimension of this algebra is `n^2`::
1157 sage: set_random_seed()
1158 sage: n = ZZ.random_element(1,5)
1159 sage: J = ComplexHermitianEJA(n)
1160 sage: J.dimension() == n^2
1163 The Jordan multiplication is what we think it is::
1165 sage: set_random_seed()
1166 sage: n = ZZ.random_element(1,5)
1167 sage: J = ComplexHermitianEJA(n)
1168 sage: x = J.random_element()
1169 sage: y = J.random_element()
1170 sage: actual = (x*y).natural_representation()
1171 sage: X = x.natural_representation()
1172 sage: Y = y.natural_representation()
1173 sage: expected = (X*Y + Y*X)/2
1174 sage: actual == expected
1176 sage: J(expected) == x*y
1180 def __init__(self
, n
, field
=QQ
):
1181 S
= _complex_hermitian_basis(n
)
1182 Qs
= _multiplication_table_from_matrix_basis(S
)
1184 fdeja
= super(ComplexHermitianEJA
, self
)
1185 return fdeja
.__init
__(field
,
1191 def inner_product(self
, x
, y
):
1192 # Since a+bi on the diagonal is represented as
1197 # we'll double-count the "a" entries if we take the trace of
1199 return _matrix_ip(x
,y
)/2
1202 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1204 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1205 matrices, the usual symmetric Jordan product, and the
1206 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1211 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1215 The dimension of this algebra is `n^2`::
1217 sage: set_random_seed()
1218 sage: n = ZZ.random_element(1,5)
1219 sage: J = QuaternionHermitianEJA(n)
1220 sage: J.dimension() == 2*(n^2) - n
1223 The Jordan multiplication is what we think it is::
1225 sage: set_random_seed()
1226 sage: n = ZZ.random_element(1,5)
1227 sage: J = QuaternionHermitianEJA(n)
1228 sage: x = J.random_element()
1229 sage: y = J.random_element()
1230 sage: actual = (x*y).natural_representation()
1231 sage: X = x.natural_representation()
1232 sage: Y = y.natural_representation()
1233 sage: expected = (X*Y + Y*X)/2
1234 sage: actual == expected
1236 sage: J(expected) == x*y
1240 def __init__(self
, n
, field
=QQ
):
1241 S
= _quaternion_hermitian_basis(n
)
1242 Qs
= _multiplication_table_from_matrix_basis(S
)
1244 fdeja
= super(QuaternionHermitianEJA
, self
)
1245 return fdeja
.__init
__(field
,
1250 def inner_product(self
, x
, y
):
1251 # Since a+bi+cj+dk on the diagonal is represented as
1253 # a + bi +cj + dk = [ a b c d]
1258 # we'll quadruple-count the "a" entries if we take the trace of
1260 return _matrix_ip(x
,y
)/4
1263 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1265 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1266 with the usual inner product and jordan product ``x*y =
1267 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1272 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1276 This multiplication table can be verified by hand::
1278 sage: J = JordanSpinEJA(4)
1279 sage: e0,e1,e2,e3 = J.gens()
1296 def __init__(self
, n
, field
=QQ
):
1297 V
= VectorSpace(field
, n
)
1298 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1308 z0
= x
.inner_product(y
)
1309 zbar
= y0
*xbar
+ x0
*ybar
1310 z
= V([z0
] + zbar
.list())
1311 mult_table
[i
][j
] = z
1313 # The rank of the spin algebra is two, unless we're in a
1314 # one-dimensional ambient space (because the rank is bounded by
1315 # the ambient dimension).
1316 fdeja
= super(JordanSpinEJA
, self
)
1317 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2))
1319 def inner_product(self
, x
, y
):
1320 return _usual_ip(x
,y
)