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
51 self
._multiplication
_table
= mult_table
53 category
= FiniteDimensionalAlgebrasWithBasis(field
).Unital()
54 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
56 range(len(mult_table
)),
59 self
.print_options(bracket
='')
62 def _element_constructor_(self
, elt
):
64 Construct an element of this algebra from its natural
67 This gets called only after the parent element _call_ method
68 fails to find a coercion for the argument.
72 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
73 ....: RealCartesianProductEJA,
74 ....: RealSymmetricEJA)
78 The identity in `S^n` is converted to the identity in the EJA::
80 sage: J = RealSymmetricEJA(3)
81 sage: I = matrix.identity(QQ,3)
85 This skew-symmetric matrix can't be represented in the EJA::
87 sage: J = RealSymmetricEJA(3)
88 sage: A = matrix(QQ,3, lambda i,j: i-j)
90 Traceback (most recent call last):
92 ArithmeticError: vector is not in free module
96 Ensure that we can convert any element of the two non-matrix
97 simple algebras (whose natural representations are their usual
98 vector representations) back and forth faithfully::
100 sage: set_random_seed()
101 sage: J = RealCartesianProductEJA(5)
102 sage: x = J.random_element()
103 sage: J(x.to_vector().column()) == x
105 sage: J = JordanSpinEJA(5)
106 sage: x = J.random_element()
107 sage: J(x.to_vector().column()) == x
111 natural_basis
= self
.natural_basis()
112 if elt
not in natural_basis
[0].matrix_space():
113 raise ValueError("not a naturally-represented algebra element")
115 # Thanks for nothing! Matrix spaces aren't vector
116 # spaces in Sage, so we have to figure out its
117 # natural-basis coordinates ourselves.
118 V
= VectorSpace(elt
.base_ring(), elt
.nrows()*elt
.ncols())
119 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
120 coords
= W
.coordinate_vector(_mat2vec(elt
))
121 return self
.from_vector(coords
)
126 Return a string representation of ``self``.
130 sage: from mjo.eja.eja_algebra import JordanSpinEJA
134 Ensure that it says what we think it says::
136 sage: JordanSpinEJA(2, field=QQ)
137 Euclidean Jordan algebra of degree 2 over Rational Field
138 sage: JordanSpinEJA(3, field=RDF)
139 Euclidean Jordan algebra of degree 3 over Real Double Field
142 # TODO: change this to say "dimension" and fix all the tests.
143 fmt
= "Euclidean Jordan algebra of degree {} over {}"
144 return fmt
.format(self
.dimension(), self
.base_ring())
146 def product_on_basis(self
, i
, j
):
149 Lei
= self
._multiplication
_table
[i
]
150 return self
.from_vector(Lei
*ej
.to_vector())
152 def _a_regular_element(self
):
154 Guess a regular element. Needed to compute the basis for our
155 characteristic polynomial coefficients.
159 sage: from mjo.eja.eja_algebra import random_eja
163 Ensure that this hacky method succeeds for every algebra that we
164 know how to construct::
166 sage: set_random_seed()
167 sage: J = random_eja()
168 sage: J._a_regular_element().is_regular()
173 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
174 if not z
.is_regular():
175 raise ValueError("don't know a regular element")
180 def _charpoly_basis_space(self
):
182 Return the vector space spanned by the basis used in our
183 characteristic polynomial coefficients. This is used not only to
184 compute those coefficients, but also any time we need to
185 evaluate the coefficients (like when we compute the trace or
188 z
= self
._a
_regular
_element
()
189 V
= self
.vector_space()
190 V1
= V
.span_of_basis( (z
**k
).to_vector() for k
in range(self
.rank()) )
191 b
= (V1
.basis() + V1
.complement().basis())
192 return V
.span_of_basis(b
)
196 def _charpoly_coeff(self
, i
):
198 Return the coefficient polynomial "a_{i}" of this algebra's
199 general characteristic polynomial.
201 Having this be a separate cached method lets us compute and
202 store the trace/determinant (a_{r-1} and a_{0} respectively)
203 separate from the entire characteristic polynomial.
205 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
206 R
= A_of_x
.base_ring()
208 # Guaranteed by theory
211 # Danger: the in-place modification is done for performance
212 # reasons (reconstructing a matrix with huge polynomial
213 # entries is slow), but I don't know how cached_method works,
214 # so it's highly possible that we're modifying some global
215 # list variable by reference, here. In other words, you
216 # probably shouldn't call this method twice on the same
217 # algebra, at the same time, in two threads
218 Ai_orig
= A_of_x
.column(i
)
219 A_of_x
.set_column(i
,xr
)
220 numerator
= A_of_x
.det()
221 A_of_x
.set_column(i
,Ai_orig
)
223 # We're relying on the theory here to ensure that each a_i is
224 # indeed back in R, and the added negative signs are to make
225 # the whole charpoly expression sum to zero.
226 return R(-numerator
/detA
)
230 def _charpoly_matrix_system(self
):
232 Compute the matrix whose entries A_ij are polynomials in
233 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
234 corresponding to `x^r` and the determinent of the matrix A =
235 [A_ij]. In other words, all of the fixed (cachable) data needed
236 to compute the coefficients of the characteristic polynomial.
241 # Construct a new algebra over a multivariate polynomial ring...
242 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
243 R
= PolynomialRing(self
.base_ring(), names
)
244 J
= FiniteDimensionalEuclideanJordanAlgebra(
246 tuple(self
._multiplication
_table
),
249 idmat
= matrix
.identity(J
.base_ring(), n
)
251 W
= self
._charpoly
_basis
_space
()
252 W
= W
.change_ring(R
.fraction_field())
254 # Starting with the standard coordinates x = (X1,X2,...,Xn)
255 # and then converting the entries to W-coordinates allows us
256 # to pass in the standard coordinates to the charpoly and get
257 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
260 # W.coordinates(x^2) eval'd at (standard z-coords)
264 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
266 # We want the middle equivalent thing in our matrix, but use
267 # the first equivalent thing instead so that we can pass in
268 # standard coordinates.
269 x
= J
.from_vector(W(R
.gens()))
271 # Handle the zeroth power separately, because computing
272 # the unit element in J is mathematically suspect.
273 x0
= W
.coordinate_vector(self
.one().to_vector())
275 l1
+= [ W
.coordinate_vector((x
**k
).to_vector()).column()
276 for k
in range(1,r
) ]
277 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
278 A_of_x
= matrix
.block(R
, 1, n
, (l1
+ l2
))
279 xr
= W
.coordinate_vector((x
**r
).to_vector())
280 return (A_of_x
, x
, xr
, A_of_x
.det())
284 def characteristic_polynomial(self
):
286 Return a characteristic polynomial that works for all elements
289 The resulting polynomial has `n+1` variables, where `n` is the
290 dimension of this algebra. The first `n` variables correspond to
291 the coordinates of an algebra element: when evaluated at the
292 coordinates of an algebra element with respect to a certain
293 basis, the result is a univariate polynomial (in the one
294 remaining variable ``t``), namely the characteristic polynomial
299 sage: from mjo.eja.eja_algebra import JordanSpinEJA
303 The characteristic polynomial in the spin algebra is given in
304 Alizadeh, Example 11.11::
306 sage: J = JordanSpinEJA(3)
307 sage: p = J.characteristic_polynomial(); p
308 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
309 sage: xvec = J.one().to_vector()
317 # The list of coefficient polynomials a_1, a_2, ..., a_n.
318 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
320 # We go to a bit of trouble here to reorder the
321 # indeterminates, so that it's easier to evaluate the
322 # characteristic polynomial at x's coordinates and get back
323 # something in terms of t, which is what we want.
325 S
= PolynomialRing(self
.base_ring(),'t')
327 S
= PolynomialRing(S
, R
.variable_names())
330 # Note: all entries past the rth should be zero. The
331 # coefficient of the highest power (x^r) is 1, but it doesn't
332 # appear in the solution vector which contains coefficients
333 # for the other powers (to make them sum to x^r).
335 a
[r
] = 1 # corresponds to x^r
337 # When the rank is equal to the dimension, trying to
338 # assign a[r] goes out-of-bounds.
339 a
.append(1) # corresponds to x^r
341 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
344 def inner_product(self
, x
, y
):
346 The inner product associated with this Euclidean Jordan algebra.
348 Defaults to the trace inner product, but can be overridden by
349 subclasses if they are sure that the necessary properties are
354 sage: from mjo.eja.eja_algebra import random_eja
358 The inner product must satisfy its axiom for this algebra to truly
359 be a Euclidean Jordan Algebra::
361 sage: set_random_seed()
362 sage: J = random_eja()
363 sage: x = J.random_element()
364 sage: y = J.random_element()
365 sage: z = J.random_element()
366 sage: (x*y).inner_product(z) == y.inner_product(x*z)
370 if (not x
in self
) or (not y
in self
):
371 raise TypeError("arguments must live in this algebra")
372 return x
.trace_inner_product(y
)
375 def natural_basis(self
):
377 Return a more-natural representation of this algebra's basis.
379 Every finite-dimensional Euclidean Jordan Algebra is a direct
380 sum of five simple algebras, four of which comprise Hermitian
381 matrices. This method returns the original "natural" basis
382 for our underlying vector space. (Typically, the natural basis
383 is used to construct the multiplication table in the first place.)
385 Note that this will always return a matrix. The standard basis
386 in `R^n` will be returned as `n`-by-`1` column matrices.
390 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
391 ....: RealSymmetricEJA)
395 sage: J = RealSymmetricEJA(2)
397 Finite family {0: e0, 1: e1, 2: e2}
398 sage: J.natural_basis()
406 sage: J = JordanSpinEJA(2)
408 Finite family {0: e0, 1: e1}
409 sage: J.natural_basis()
416 if self
._natural
_basis
is None:
417 return tuple( b
.to_vector().column() for b
in self
.basis() )
419 return self
._natural
_basis
425 Return the unit element of this algebra.
429 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
434 sage: J = RealCartesianProductEJA(5)
436 e0 + e1 + e2 + e3 + e4
440 The identity element acts like the identity::
442 sage: set_random_seed()
443 sage: J = random_eja()
444 sage: x = J.random_element()
445 sage: J.one()*x == x and x*J.one() == x
448 The matrix of the unit element's operator is the identity::
450 sage: set_random_seed()
451 sage: J = random_eja()
452 sage: actual = J.one().operator().matrix()
453 sage: expected = matrix.identity(J.base_ring(), J.dimension())
454 sage: actual == expected
458 # We can brute-force compute the matrices of the operators
459 # that correspond to the basis elements of this algebra.
460 # If some linear combination of those basis elements is the
461 # algebra identity, then the same linear combination of
462 # their matrices has to be the identity matrix.
464 # Of course, matrices aren't vectors in sage, so we have to
465 # appeal to the "long vectors" isometry.
466 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
468 # Now we use basis linear algebra to find the coefficients,
469 # of the matrices-as-vectors-linear-combination, which should
470 # work for the original algebra basis too.
471 A
= matrix
.column(self
.base_ring(), oper_vecs
)
473 # We used the isometry on the left-hand side already, but we
474 # still need to do it for the right-hand side. Recall that we
475 # wanted something that summed to the identity matrix.
476 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
478 # Now if there's an identity element in the algebra, this should work.
479 coeffs
= A
.solve_right(b
)
480 return self
.linear_combination(zip(self
.gens(), coeffs
))
485 Return the rank of this EJA.
489 The author knows of no algorithm to compute the rank of an EJA
490 where only the multiplication table is known. In lieu of one, we
491 require the rank to be specified when the algebra is created,
492 and simply pass along that number here.
496 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
497 ....: RealSymmetricEJA,
498 ....: ComplexHermitianEJA,
499 ....: QuaternionHermitianEJA,
504 The rank of the Jordan spin algebra is always two::
506 sage: JordanSpinEJA(2).rank()
508 sage: JordanSpinEJA(3).rank()
510 sage: JordanSpinEJA(4).rank()
513 The rank of the `n`-by-`n` Hermitian real, complex, or
514 quaternion matrices is `n`::
516 sage: RealSymmetricEJA(2).rank()
518 sage: ComplexHermitianEJA(2).rank()
520 sage: QuaternionHermitianEJA(2).rank()
522 sage: RealSymmetricEJA(5).rank()
524 sage: ComplexHermitianEJA(5).rank()
526 sage: QuaternionHermitianEJA(5).rank()
531 Ensure that every EJA that we know how to construct has a
532 positive integer rank::
534 sage: set_random_seed()
535 sage: r = random_eja().rank()
536 sage: r in ZZ and r > 0
543 def vector_space(self
):
545 Return the vector space that underlies this algebra.
549 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
553 sage: J = RealSymmetricEJA(2)
554 sage: J.vector_space()
555 Vector space of dimension 3 over Rational Field
558 return self
.zero().to_vector().parent().ambient_vector_space()
561 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
564 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
566 Return the Euclidean Jordan Algebra corresponding to the set
567 `R^n` under the Hadamard product.
569 Note: this is nothing more than the Cartesian product of ``n``
570 copies of the spin algebra. Once Cartesian product algebras
571 are implemented, this can go.
575 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
579 This multiplication table can be verified by hand::
581 sage: J = RealCartesianProductEJA(3)
582 sage: e0,e1,e2 = J.gens()
597 def __init__(self
, n
, field
=QQ
):
598 # The superclass constructor takes a list of matrices, the ith
599 # representing right multiplication by the ith basis element
600 # in the vector space. So if e_1 = (1,0,0), then right
601 # (Hadamard) multiplication of x by e_1 picks out the first
602 # component of x; and likewise for the ith basis element e_i.
603 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
606 fdeja
= super(RealCartesianProductEJA
, self
)
607 return fdeja
.__init
__(field
, Qs
, rank
=n
)
609 def inner_product(self
, x
, y
):
610 return _usual_ip(x
,y
)
615 Return a "random" finite-dimensional Euclidean Jordan Algebra.
619 For now, we choose a random natural number ``n`` (greater than zero)
620 and then give you back one of the following:
622 * The cartesian product of the rational numbers ``n`` times; this is
623 ``QQ^n`` with the Hadamard product.
625 * The Jordan spin algebra on ``QQ^n``.
627 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
630 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
631 in the space of ``2n``-by-``2n`` real symmetric matrices.
633 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
634 in the space of ``4n``-by-``4n`` real symmetric matrices.
636 Later this might be extended to return Cartesian products of the
641 sage: from mjo.eja.eja_algebra import random_eja
646 Euclidean Jordan algebra of degree...
650 # The max_n component lets us choose different upper bounds on the
651 # value "n" that gets passed to the constructor. This is needed
652 # because e.g. R^{10} is reasonable to test, while the Hermitian
653 # 10-by-10 quaternion matrices are not.
654 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
656 (RealSymmetricEJA
, 5),
657 (ComplexHermitianEJA
, 4),
658 (QuaternionHermitianEJA
, 3)])
659 n
= ZZ
.random_element(1, max_n
)
660 return constructor(n
, field
=QQ
)
664 def _real_symmetric_basis(n
, field
=QQ
):
666 Return a basis for the space of real symmetric n-by-n matrices.
668 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
672 for j
in xrange(i
+1):
673 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
677 # Beware, orthogonal but not normalized!
678 Sij
= Eij
+ Eij
.transpose()
683 def _complex_hermitian_basis(n
, field
=QQ
):
685 Returns a basis for the space of complex Hermitian n-by-n matrices.
689 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
693 sage: set_random_seed()
694 sage: n = ZZ.random_element(1,5)
695 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
699 F
= QuadraticField(-1, 'I')
702 # This is like the symmetric case, but we need to be careful:
704 # * We want conjugate-symmetry, not just symmetry.
705 # * The diagonal will (as a result) be real.
709 for j
in xrange(i
+1):
710 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
712 Sij
= _embed_complex_matrix(Eij
)
715 # Beware, orthogonal but not normalized! The second one
716 # has a minus because it's conjugated.
717 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
719 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
724 def _quaternion_hermitian_basis(n
, field
=QQ
):
726 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
730 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
734 sage: set_random_seed()
735 sage: n = ZZ.random_element(1,5)
736 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
740 Q
= QuaternionAlgebra(QQ
,-1,-1)
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(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
753 Sij
= _embed_quaternion_matrix(Eij
)
756 # Beware, orthogonal but not normalized! The second,
757 # third, and fourth ones have a minus because they're
759 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
761 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
763 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
765 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
771 def _multiplication_table_from_matrix_basis(basis
):
773 At least three of the five simple Euclidean Jordan algebras have the
774 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
775 multiplication on the right is matrix multiplication. Given a basis
776 for the underlying matrix space, this function returns a
777 multiplication table (obtained by looping through the basis
778 elements) for an algebra of those matrices. A reordered copy
779 of the basis is also returned to work around the fact that
780 the ``span()`` in this function will change the order of the basis
781 from what we think it is, to... something else.
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
)
796 # Brute force the multiplication-by-s matrix by looping
797 # through all elements of the basis and doing the computation
798 # to find out what the corresponding row should be.
801 this_col
= _mat2vec((s
*t
+ t
*s
)/2)
802 Q_cols
.append(W
.coordinates(this_col
))
803 Q
= matrix
.column(field
, W
.dimension(), Q_cols
)
809 def _embed_complex_matrix(M
):
811 Embed the n-by-n complex matrix ``M`` into the space of real
812 matrices of size 2n-by-2n via the map the sends each entry `z = a +
813 bi` to the block matrix ``[[a,b],[-b,a]]``.
817 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
821 sage: F = QuadraticField(-1,'i')
822 sage: x1 = F(4 - 2*i)
823 sage: x2 = F(1 + 2*i)
826 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
827 sage: _embed_complex_matrix(M)
836 Embedding is a homomorphism (isomorphism, in fact)::
838 sage: set_random_seed()
839 sage: n = ZZ.random_element(5)
840 sage: F = QuadraticField(-1, 'i')
841 sage: X = random_matrix(F, n)
842 sage: Y = random_matrix(F, n)
843 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
844 sage: expected = _embed_complex_matrix(X*Y)
845 sage: actual == expected
851 raise ValueError("the matrix 'M' must be square")
852 field
= M
.base_ring()
857 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
859 # We can drop the imaginaries here.
860 return matrix
.block(field
.base_ring(), n
, blocks
)
863 def _unembed_complex_matrix(M
):
865 The inverse of _embed_complex_matrix().
869 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
870 ....: _unembed_complex_matrix)
874 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
875 ....: [-2, 1, -4, 3],
876 ....: [ 9, 10, 11, 12],
877 ....: [-10, 9, -12, 11] ])
878 sage: _unembed_complex_matrix(A)
880 [ 10*i + 9 12*i + 11]
884 Unembedding is the inverse of embedding::
886 sage: set_random_seed()
887 sage: F = QuadraticField(-1, 'i')
888 sage: M = random_matrix(F, 3)
889 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
895 raise ValueError("the matrix 'M' must be square")
896 if not n
.mod(2).is_zero():
897 raise ValueError("the matrix 'M' must be a complex embedding")
899 F
= QuadraticField(-1, 'i')
902 # Go top-left to bottom-right (reading order), converting every
903 # 2-by-2 block we see to a single complex element.
905 for k
in xrange(n
/2):
906 for j
in xrange(n
/2):
907 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
908 if submat
[0,0] != submat
[1,1]:
909 raise ValueError('bad on-diagonal submatrix')
910 if submat
[0,1] != -submat
[1,0]:
911 raise ValueError('bad off-diagonal submatrix')
912 z
= submat
[0,0] + submat
[0,1]*i
915 return matrix(F
, n
/2, elements
)
918 def _embed_quaternion_matrix(M
):
920 Embed the n-by-n quaternion matrix ``M`` into the space of real
921 matrices of size 4n-by-4n by first sending each quaternion entry
922 `z = a + bi + cj + dk` to the block-complex matrix
923 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
928 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
932 sage: Q = QuaternionAlgebra(QQ,-1,-1)
933 sage: i,j,k = Q.gens()
934 sage: x = 1 + 2*i + 3*j + 4*k
935 sage: M = matrix(Q, 1, [[x]])
936 sage: _embed_quaternion_matrix(M)
942 Embedding is a homomorphism (isomorphism, in fact)::
944 sage: set_random_seed()
945 sage: n = ZZ.random_element(5)
946 sage: Q = QuaternionAlgebra(QQ,-1,-1)
947 sage: X = random_matrix(Q, n)
948 sage: Y = random_matrix(Q, n)
949 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
950 sage: expected = _embed_quaternion_matrix(X*Y)
951 sage: actual == expected
955 quaternions
= M
.base_ring()
958 raise ValueError("the matrix 'M' must be square")
960 F
= QuadraticField(-1, 'i')
965 t
= z
.coefficient_tuple()
970 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
971 [-c
+ d
*i
, a
- b
*i
]])
972 blocks
.append(_embed_complex_matrix(cplx_matrix
))
974 # We should have real entries by now, so use the realest field
975 # we've got for the return value.
976 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
979 def _unembed_quaternion_matrix(M
):
981 The inverse of _embed_quaternion_matrix().
985 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
986 ....: _unembed_quaternion_matrix)
990 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
991 ....: [-2, 1, -4, 3],
992 ....: [-3, 4, 1, -2],
993 ....: [-4, -3, 2, 1]])
994 sage: _unembed_quaternion_matrix(M)
995 [1 + 2*i + 3*j + 4*k]
999 Unembedding is the inverse of embedding::
1001 sage: set_random_seed()
1002 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1003 sage: M = random_matrix(Q, 3)
1004 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1010 raise ValueError("the matrix 'M' must be square")
1011 if not n
.mod(4).is_zero():
1012 raise ValueError("the matrix 'M' must be a complex embedding")
1014 Q
= QuaternionAlgebra(QQ
,-1,-1)
1017 # Go top-left to bottom-right (reading order), converting every
1018 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1021 for l
in xrange(n
/4):
1022 for m
in xrange(n
/4):
1023 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1024 if submat
[0,0] != submat
[1,1].conjugate():
1025 raise ValueError('bad on-diagonal submatrix')
1026 if submat
[0,1] != -submat
[1,0].conjugate():
1027 raise ValueError('bad off-diagonal submatrix')
1028 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1029 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1032 return matrix(Q
, n
/4, elements
)
1035 # The usual inner product on R^n.
1037 return x
.to_vector().inner_product(y
.to_vector())
1039 # The inner product used for the real symmetric simple EJA.
1040 # We keep it as a separate function because e.g. the complex
1041 # algebra uses the same inner product, except divided by 2.
1042 def _matrix_ip(X
,Y
):
1043 X_mat
= X
.natural_representation()
1044 Y_mat
= Y
.natural_representation()
1045 return (X_mat
*Y_mat
).trace()
1048 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1050 The rank-n simple EJA consisting of real symmetric n-by-n
1051 matrices, the usual symmetric Jordan product, and the trace inner
1052 product. It has dimension `(n^2 + n)/2` over the reals.
1056 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1060 sage: J = RealSymmetricEJA(2)
1061 sage: e0, e1, e2 = J.gens()
1071 The dimension of this algebra is `(n^2 + n) / 2`::
1073 sage: set_random_seed()
1074 sage: n = ZZ.random_element(1,5)
1075 sage: J = RealSymmetricEJA(n)
1076 sage: J.dimension() == (n^2 + n)/2
1079 The Jordan multiplication is what we think it is::
1081 sage: set_random_seed()
1082 sage: n = ZZ.random_element(1,5)
1083 sage: J = RealSymmetricEJA(n)
1084 sage: x = J.random_element()
1085 sage: y = J.random_element()
1086 sage: actual = (x*y).natural_representation()
1087 sage: X = x.natural_representation()
1088 sage: Y = y.natural_representation()
1089 sage: expected = (X*Y + Y*X)/2
1090 sage: actual == expected
1092 sage: J(expected) == x*y
1096 def __init__(self
, n
, field
=QQ
):
1097 S
= _real_symmetric_basis(n
, field
=field
)
1098 Qs
= _multiplication_table_from_matrix_basis(S
)
1100 fdeja
= super(RealSymmetricEJA
, self
)
1101 return fdeja
.__init
__(field
,
1106 def inner_product(self
, x
, y
):
1107 return _matrix_ip(x
,y
)
1110 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1112 The rank-n simple EJA consisting of complex Hermitian n-by-n
1113 matrices over the real numbers, the usual symmetric Jordan product,
1114 and the real-part-of-trace inner product. It has dimension `n^2` over
1119 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1123 The dimension of this algebra is `n^2`::
1125 sage: set_random_seed()
1126 sage: n = ZZ.random_element(1,5)
1127 sage: J = ComplexHermitianEJA(n)
1128 sage: J.dimension() == n^2
1131 The Jordan multiplication is what we think it is::
1133 sage: set_random_seed()
1134 sage: n = ZZ.random_element(1,5)
1135 sage: J = ComplexHermitianEJA(n)
1136 sage: x = J.random_element()
1137 sage: y = J.random_element()
1138 sage: actual = (x*y).natural_representation()
1139 sage: X = x.natural_representation()
1140 sage: Y = y.natural_representation()
1141 sage: expected = (X*Y + Y*X)/2
1142 sage: actual == expected
1144 sage: J(expected) == x*y
1148 def __init__(self
, n
, field
=QQ
):
1149 S
= _complex_hermitian_basis(n
)
1150 Qs
= _multiplication_table_from_matrix_basis(S
)
1152 fdeja
= super(ComplexHermitianEJA
, self
)
1153 return fdeja
.__init
__(field
,
1159 def inner_product(self
, x
, y
):
1160 # Since a+bi on the diagonal is represented as
1165 # we'll double-count the "a" entries if we take the trace of
1167 return _matrix_ip(x
,y
)/2
1170 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1172 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1173 matrices, the usual symmetric Jordan product, and the
1174 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1179 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1183 The dimension of this algebra is `n^2`::
1185 sage: set_random_seed()
1186 sage: n = ZZ.random_element(1,5)
1187 sage: J = QuaternionHermitianEJA(n)
1188 sage: J.dimension() == 2*(n^2) - n
1191 The Jordan multiplication is what we think it is::
1193 sage: set_random_seed()
1194 sage: n = ZZ.random_element(1,5)
1195 sage: J = QuaternionHermitianEJA(n)
1196 sage: x = J.random_element()
1197 sage: y = J.random_element()
1198 sage: actual = (x*y).natural_representation()
1199 sage: X = x.natural_representation()
1200 sage: Y = y.natural_representation()
1201 sage: expected = (X*Y + Y*X)/2
1202 sage: actual == expected
1204 sage: J(expected) == x*y
1208 def __init__(self
, n
, field
=QQ
):
1209 S
= _quaternion_hermitian_basis(n
)
1210 Qs
= _multiplication_table_from_matrix_basis(S
)
1212 fdeja
= super(QuaternionHermitianEJA
, self
)
1213 return fdeja
.__init
__(field
,
1218 def inner_product(self
, x
, y
):
1219 # Since a+bi+cj+dk on the diagonal is represented as
1221 # a + bi +cj + dk = [ a b c d]
1226 # we'll quadruple-count the "a" entries if we take the trace of
1228 return _matrix_ip(x
,y
)/4
1231 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1233 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1234 with the usual inner product and jordan product ``x*y =
1235 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1240 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1244 This multiplication table can be verified by hand::
1246 sage: J = JordanSpinEJA(4)
1247 sage: e0,e1,e2,e3 = J.gens()
1264 def __init__(self
, n
, field
=QQ
):
1266 id_matrix
= matrix
.identity(field
, n
)
1268 ei
= id_matrix
.column(i
)
1269 Qi
= matrix
.zero(field
, n
)
1271 Qi
.set_column(0, ei
)
1272 Qi
+= matrix
.diagonal(n
, [ei
[0]]*n
)
1273 # The addition of the diagonal matrix adds an extra ei[0] in the
1274 # upper-left corner of the matrix.
1275 Qi
[0,0] = Qi
[0,0] * ~
field(2)
1278 # The rank of the spin algebra is two, unless we're in a
1279 # one-dimensional ambient space (because the rank is bounded by
1280 # the ambient dimension).
1281 fdeja
= super(JordanSpinEJA
, self
)
1282 return fdeja
.__init
__(field
, Qs
, rank
=min(n
,2))
1284 def inner_product(self
, x
, y
):
1285 return _usual_ip(x
,y
)