2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 #from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra
9 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
10 from sage
.categories
.finite_dimensional_algebras_with_basis
import FiniteDimensionalAlgebrasWithBasis
11 from sage
.combinat
.free_module
import CombinatorialFreeModule
12 from sage
.matrix
.constructor
import matrix
13 from sage
.misc
.cachefunc
import cached_method
14 from sage
.misc
.prandom
import choice
15 from sage
.modules
.free_module
import VectorSpace
16 from sage
.rings
.integer_ring
import ZZ
17 from sage
.rings
.number_field
.number_field
import QuadraticField
18 from sage
.rings
.polynomial
.polynomial_ring_constructor
import PolynomialRing
19 from sage
.rings
.rational_field
import QQ
20 from sage
.structure
.element
import is_Matrix
21 from sage
.structure
.category_object
import normalize_names
23 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
24 from mjo
.eja
.eja_utils
import _mat2vec
26 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
37 sage: from mjo.eja.eja_algebra import random_eja
41 By definition, Jordan multiplication commutes::
43 sage: set_random_seed()
44 sage: J = random_eja()
45 sage: x = J.random_element()
46 sage: y = J.random_element()
52 self
._natural
_basis
= natural_basis
53 self
._multiplication
_table
= mult_table
55 category
= FiniteDimensionalAlgebrasWithBasis(field
).Unital()
56 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
58 range(len(mult_table
)),
61 self
.print_options(bracket
='')
64 def _element_constructor_(self
, elt
):
66 Construct an element of this algebra from its natural
69 This gets called only after the parent element _call_ method
70 fails to find a coercion for the argument.
74 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
75 ....: RealCartesianProductEJA,
76 ....: RealSymmetricEJA)
80 The identity in `S^n` is converted to the identity in the EJA::
82 sage: J = RealSymmetricEJA(3)
83 sage: I = matrix.identity(QQ,3)
87 This skew-symmetric matrix can't be represented in the EJA::
89 sage: J = RealSymmetricEJA(3)
90 sage: A = matrix(QQ,3, lambda i,j: i-j)
92 Traceback (most recent call last):
94 ArithmeticError: vector is not in free module
98 Ensure that we can convert any element of the two non-matrix
99 simple algebras (whose natural representations are their usual
100 vector representations) back and forth faithfully::
102 sage: set_random_seed()
103 sage: J = RealCartesianProductEJA(5)
104 sage: x = J.random_element()
105 sage: J(x.to_vector().column()) == x
107 sage: J = JordanSpinEJA(5)
108 sage: x = J.random_element()
109 sage: J(x.to_vector().column()) == x
113 natural_basis
= self
.natural_basis()
114 if elt
not in natural_basis
[0].matrix_space():
115 raise ValueError("not a naturally-represented algebra element")
117 # Thanks for nothing! Matrix spaces aren't vector
118 # spaces in Sage, so we have to figure out its
119 # natural-basis coordinates ourselves.
120 V
= VectorSpace(elt
.base_ring(), elt
.nrows()*elt
.ncols())
121 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
122 coords
= W
.coordinate_vector(_mat2vec(elt
))
123 return self
.from_vector(coords
)
128 Return a string representation of ``self``.
132 sage: from mjo.eja.eja_algebra import JordanSpinEJA
136 Ensure that it says what we think it says::
138 sage: JordanSpinEJA(2, field=QQ)
139 Euclidean Jordan algebra of degree 2 over Rational Field
140 sage: JordanSpinEJA(3, field=RDF)
141 Euclidean Jordan algebra of degree 3 over Real Double Field
144 # TODO: change this to say "dimension" and fix all the tests.
145 fmt
= "Euclidean Jordan algebra of degree {} over {}"
146 return fmt
.format(self
.dimension(), self
.base_ring())
148 def product_on_basis(self
, i
, j
):
151 Lei
= self
._multiplication
_table
[i
]
152 return self
.from_vector(Lei
*ej
.to_vector())
154 def _a_regular_element(self
):
156 Guess a regular element. Needed to compute the basis for our
157 characteristic polynomial coefficients.
161 sage: from mjo.eja.eja_algebra import random_eja
165 Ensure that this hacky method succeeds for every algebra that we
166 know how to construct::
168 sage: set_random_seed()
169 sage: J = random_eja()
170 sage: J._a_regular_element().is_regular()
175 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
176 if not z
.is_regular():
177 raise ValueError("don't know a regular element")
182 def _charpoly_basis_space(self
):
184 Return the vector space spanned by the basis used in our
185 characteristic polynomial coefficients. This is used not only to
186 compute those coefficients, but also any time we need to
187 evaluate the coefficients (like when we compute the trace or
190 z
= self
._a
_regular
_element
()
191 V
= self
.vector_space()
192 V1
= V
.span_of_basis( (z
**k
).to_vector() for k
in range(self
.rank()) )
193 b
= (V1
.basis() + V1
.complement().basis())
194 return V
.span_of_basis(b
)
198 def _charpoly_coeff(self
, i
):
200 Return the coefficient polynomial "a_{i}" of this algebra's
201 general characteristic polynomial.
203 Having this be a separate cached method lets us compute and
204 store the trace/determinant (a_{r-1} and a_{0} respectively)
205 separate from the entire characteristic polynomial.
207 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
208 R
= A_of_x
.base_ring()
210 # Guaranteed by theory
213 # Danger: the in-place modification is done for performance
214 # reasons (reconstructing a matrix with huge polynomial
215 # entries is slow), but I don't know how cached_method works,
216 # so it's highly possible that we're modifying some global
217 # list variable by reference, here. In other words, you
218 # probably shouldn't call this method twice on the same
219 # algebra, at the same time, in two threads
220 Ai_orig
= A_of_x
.column(i
)
221 A_of_x
.set_column(i
,xr
)
222 numerator
= A_of_x
.det()
223 A_of_x
.set_column(i
,Ai_orig
)
225 # We're relying on the theory here to ensure that each a_i is
226 # indeed back in R, and the added negative signs are to make
227 # the whole charpoly expression sum to zero.
228 return R(-numerator
/detA
)
232 def _charpoly_matrix_system(self
):
234 Compute the matrix whose entries A_ij are polynomials in
235 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
236 corresponding to `x^r` and the determinent of the matrix A =
237 [A_ij]. In other words, all of the fixed (cachable) data needed
238 to compute the coefficients of the characteristic polynomial.
243 # Construct a new algebra over a multivariate polynomial ring...
244 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
245 R
= PolynomialRing(self
.base_ring(), names
)
246 J
= FiniteDimensionalEuclideanJordanAlgebra(
248 tuple(self
._multiplication
_table
),
251 idmat
= matrix
.identity(J
.base_ring(), n
)
253 W
= self
._charpoly
_basis
_space
()
254 W
= W
.change_ring(R
.fraction_field())
256 # Starting with the standard coordinates x = (X1,X2,...,Xn)
257 # and then converting the entries to W-coordinates allows us
258 # to pass in the standard coordinates to the charpoly and get
259 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
262 # W.coordinates(x^2) eval'd at (standard z-coords)
266 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
268 # We want the middle equivalent thing in our matrix, but use
269 # the first equivalent thing instead so that we can pass in
270 # standard coordinates.
271 x
= J
.from_vector(W(R
.gens()))
273 # Handle the zeroth power separately, because computing
274 # the unit element in J is mathematically suspect.
275 x0
= W
.coordinate_vector(self
.one().to_vector())
277 l1
+= [ W
.coordinate_vector((x
**k
).to_vector()).column()
278 for k
in range(1,r
) ]
279 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
280 A_of_x
= matrix
.block(R
, 1, n
, (l1
+ l2
))
281 xr
= W
.coordinate_vector((x
**r
).to_vector())
282 return (A_of_x
, x
, xr
, A_of_x
.det())
286 def characteristic_polynomial(self
):
288 Return a characteristic polynomial that works for all elements
291 The resulting polynomial has `n+1` variables, where `n` is the
292 dimension of this algebra. The first `n` variables correspond to
293 the coordinates of an algebra element: when evaluated at the
294 coordinates of an algebra element with respect to a certain
295 basis, the result is a univariate polynomial (in the one
296 remaining variable ``t``), namely the characteristic polynomial
301 sage: from mjo.eja.eja_algebra import JordanSpinEJA
305 The characteristic polynomial in the spin algebra is given in
306 Alizadeh, Example 11.11::
308 sage: J = JordanSpinEJA(3)
309 sage: p = J.characteristic_polynomial(); p
310 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
311 sage: xvec = J.one().to_vector()
319 # The list of coefficient polynomials a_1, a_2, ..., a_n.
320 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
322 # We go to a bit of trouble here to reorder the
323 # indeterminates, so that it's easier to evaluate the
324 # characteristic polynomial at x's coordinates and get back
325 # something in terms of t, which is what we want.
327 S
= PolynomialRing(self
.base_ring(),'t')
329 S
= PolynomialRing(S
, R
.variable_names())
332 # Note: all entries past the rth should be zero. The
333 # coefficient of the highest power (x^r) is 1, but it doesn't
334 # appear in the solution vector which contains coefficients
335 # for the other powers (to make them sum to x^r).
337 a
[r
] = 1 # corresponds to x^r
339 # When the rank is equal to the dimension, trying to
340 # assign a[r] goes out-of-bounds.
341 a
.append(1) # corresponds to x^r
343 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
346 def inner_product(self
, x
, y
):
348 The inner product associated with this Euclidean Jordan algebra.
350 Defaults to the trace inner product, but can be overridden by
351 subclasses if they are sure that the necessary properties are
356 sage: from mjo.eja.eja_algebra import random_eja
360 The inner product must satisfy its axiom for this algebra to truly
361 be a Euclidean Jordan Algebra::
363 sage: set_random_seed()
364 sage: J = random_eja()
365 sage: x = J.random_element()
366 sage: y = J.random_element()
367 sage: z = J.random_element()
368 sage: (x*y).inner_product(z) == y.inner_product(x*z)
372 if (not x
in self
) or (not y
in self
):
373 raise TypeError("arguments must live in this algebra")
374 return x
.trace_inner_product(y
)
377 def natural_basis(self
):
379 Return a more-natural representation of this algebra's basis.
381 Every finite-dimensional Euclidean Jordan Algebra is a direct
382 sum of five simple algebras, four of which comprise Hermitian
383 matrices. This method returns the original "natural" basis
384 for our underlying vector space. (Typically, the natural basis
385 is used to construct the multiplication table in the first place.)
387 Note that this will always return a matrix. The standard basis
388 in `R^n` will be returned as `n`-by-`1` column matrices.
392 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
393 ....: RealSymmetricEJA)
397 sage: J = RealSymmetricEJA(2)
399 Finite family {0: e0, 1: e1, 2: e2}
400 sage: J.natural_basis()
408 sage: J = JordanSpinEJA(2)
410 Finite family {0: e0, 1: e1}
411 sage: J.natural_basis()
418 if self
._natural
_basis
is None:
419 return tuple( b
.to_vector().column() for b
in self
.basis() )
421 return self
._natural
_basis
427 Return the unit element of this algebra.
431 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
436 sage: J = RealCartesianProductEJA(5)
438 e0 + e1 + e2 + e3 + e4
442 The identity element acts like the identity::
444 sage: set_random_seed()
445 sage: J = random_eja()
446 sage: x = J.random_element()
447 sage: J.one()*x == x and x*J.one() == x
450 The matrix of the unit element's operator is the identity::
452 sage: set_random_seed()
453 sage: J = random_eja()
454 sage: actual = J.one().operator().matrix()
455 sage: expected = matrix.identity(J.base_ring(), J.dimension())
456 sage: actual == expected
460 # We can brute-force compute the matrices of the operators
461 # that correspond to the basis elements of this algebra.
462 # If some linear combination of those basis elements is the
463 # algebra identity, then the same linear combination of
464 # their matrices has to be the identity matrix.
466 # Of course, matrices aren't vectors in sage, so we have to
467 # appeal to the "long vectors" isometry.
468 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
470 # Now we use basis linear algebra to find the coefficients,
471 # of the matrices-as-vectors-linear-combination, which should
472 # work for the original algebra basis too.
473 A
= matrix
.column(self
.base_ring(), oper_vecs
)
475 # We used the isometry on the left-hand side already, but we
476 # still need to do it for the right-hand side. Recall that we
477 # wanted something that summed to the identity matrix.
478 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
480 # Now if there's an identity element in the algebra, this should work.
481 coeffs
= A
.solve_right(b
)
482 return self
.linear_combination(zip(self
.gens(), coeffs
))
487 Return the rank of this EJA.
491 The author knows of no algorithm to compute the rank of an EJA
492 where only the multiplication table is known. In lieu of one, we
493 require the rank to be specified when the algebra is created,
494 and simply pass along that number here.
498 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
499 ....: RealSymmetricEJA,
500 ....: ComplexHermitianEJA,
501 ....: QuaternionHermitianEJA,
506 The rank of the Jordan spin algebra is always two::
508 sage: JordanSpinEJA(2).rank()
510 sage: JordanSpinEJA(3).rank()
512 sage: JordanSpinEJA(4).rank()
515 The rank of the `n`-by-`n` Hermitian real, complex, or
516 quaternion matrices is `n`::
518 sage: RealSymmetricEJA(2).rank()
520 sage: ComplexHermitianEJA(2).rank()
522 sage: QuaternionHermitianEJA(2).rank()
524 sage: RealSymmetricEJA(5).rank()
526 sage: ComplexHermitianEJA(5).rank()
528 sage: QuaternionHermitianEJA(5).rank()
533 Ensure that every EJA that we know how to construct has a
534 positive integer rank::
536 sage: set_random_seed()
537 sage: r = random_eja().rank()
538 sage: r in ZZ and r > 0
545 def vector_space(self
):
547 Return the vector space that underlies this algebra.
551 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
555 sage: J = RealSymmetricEJA(2)
556 sage: J.vector_space()
557 Vector space of dimension 3 over Rational Field
560 return self
.zero().to_vector().parent().ambient_vector_space()
563 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
566 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
568 Return the Euclidean Jordan Algebra corresponding to the set
569 `R^n` under the Hadamard product.
571 Note: this is nothing more than the Cartesian product of ``n``
572 copies of the spin algebra. Once Cartesian product algebras
573 are implemented, this can go.
577 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
581 This multiplication table can be verified by hand::
583 sage: J = RealCartesianProductEJA(3)
584 sage: e0,e1,e2 = J.gens()
599 def __init__(self
, n
, field
=QQ
):
600 # The superclass constructor takes a list of matrices, the ith
601 # representing right multiplication by the ith basis element
602 # in the vector space. So if e_1 = (1,0,0), then right
603 # (Hadamard) multiplication of x by e_1 picks out the first
604 # component of x; and likewise for the ith basis element e_i.
605 Qs
= [ matrix(field
, n
, n
, lambda k
,j
: 1*(k
== j
== i
))
608 fdeja
= super(RealCartesianProductEJA
, self
)
609 return fdeja
.__init
__(field
, Qs
, rank
=n
)
611 def inner_product(self
, x
, y
):
612 return _usual_ip(x
,y
)
617 Return a "random" finite-dimensional Euclidean Jordan Algebra.
621 For now, we choose a random natural number ``n`` (greater than zero)
622 and then give you back one of the following:
624 * The cartesian product of the rational numbers ``n`` times; this is
625 ``QQ^n`` with the Hadamard product.
627 * The Jordan spin algebra on ``QQ^n``.
629 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
632 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
633 in the space of ``2n``-by-``2n`` real symmetric matrices.
635 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
636 in the space of ``4n``-by-``4n`` real symmetric matrices.
638 Later this might be extended to return Cartesian products of the
643 sage: from mjo.eja.eja_algebra import random_eja
648 Euclidean Jordan algebra of degree...
652 # The max_n component lets us choose different upper bounds on the
653 # value "n" that gets passed to the constructor. This is needed
654 # because e.g. R^{10} is reasonable to test, while the Hermitian
655 # 10-by-10 quaternion matrices are not.
656 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
658 (RealSymmetricEJA
, 5),
659 (ComplexHermitianEJA
, 4),
660 (QuaternionHermitianEJA
, 3)])
661 n
= ZZ
.random_element(1, max_n
)
662 return constructor(n
, field
=QQ
)
666 def _real_symmetric_basis(n
, field
=QQ
):
668 Return a basis for the space of real symmetric n-by-n matrices.
670 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
674 for j
in xrange(i
+1):
675 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
679 # Beware, orthogonal but not normalized!
680 Sij
= Eij
+ Eij
.transpose()
685 def _complex_hermitian_basis(n
, field
=QQ
):
687 Returns a basis for the space of complex Hermitian n-by-n matrices.
691 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
695 sage: set_random_seed()
696 sage: n = ZZ.random_element(1,5)
697 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
701 F
= QuadraticField(-1, 'I')
704 # This is like the symmetric case, but we need to be careful:
706 # * We want conjugate-symmetry, not just symmetry.
707 # * The diagonal will (as a result) be real.
711 for j
in xrange(i
+1):
712 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
714 Sij
= _embed_complex_matrix(Eij
)
717 # Beware, orthogonal but not normalized! The second one
718 # has a minus because it's conjugated.
719 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
721 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
726 def _quaternion_hermitian_basis(n
, field
=QQ
):
728 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
732 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
736 sage: set_random_seed()
737 sage: n = ZZ.random_element(1,5)
738 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
742 Q
= QuaternionAlgebra(QQ
,-1,-1)
745 # This is like the symmetric case, but we need to be careful:
747 # * We want conjugate-symmetry, not just symmetry.
748 # * The diagonal will (as a result) be real.
752 for j
in xrange(i
+1):
753 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
755 Sij
= _embed_quaternion_matrix(Eij
)
758 # Beware, orthogonal but not normalized! The second,
759 # third, and fourth ones have a minus because they're
761 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
763 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
765 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
767 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
773 def _multiplication_table_from_matrix_basis(basis
):
775 At least three of the five simple Euclidean Jordan algebras have the
776 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
777 multiplication on the right is matrix multiplication. Given a basis
778 for the underlying matrix space, this function returns a
779 multiplication table (obtained by looping through the basis
780 elements) for an algebra of those matrices. A reordered copy
781 of the basis is also returned to work around the fact that
782 the ``span()`` in this function will change the order of the basis
783 from what we think it is, to... something else.
785 # In S^2, for example, we nominally have four coordinates even
786 # though the space is of dimension three only. The vector space V
787 # is supposed to hold the entire long vector, and the subspace W
788 # of V will be spanned by the vectors that arise from symmetric
789 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
790 field
= basis
[0].base_ring()
791 dimension
= basis
[0].nrows()
793 V
= VectorSpace(field
, dimension
**2)
794 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
798 # Brute force the multiplication-by-s matrix by looping
799 # through all elements of the basis and doing the computation
800 # to find out what the corresponding row should be.
803 this_col
= _mat2vec((s
*t
+ t
*s
)/2)
804 Q_cols
.append(W
.coordinates(this_col
))
805 Q
= matrix
.column(field
, W
.dimension(), Q_cols
)
811 def _embed_complex_matrix(M
):
813 Embed the n-by-n complex matrix ``M`` into the space of real
814 matrices of size 2n-by-2n via the map the sends each entry `z = a +
815 bi` to the block matrix ``[[a,b],[-b,a]]``.
819 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
823 sage: F = QuadraticField(-1,'i')
824 sage: x1 = F(4 - 2*i)
825 sage: x2 = F(1 + 2*i)
828 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
829 sage: _embed_complex_matrix(M)
838 Embedding is a homomorphism (isomorphism, in fact)::
840 sage: set_random_seed()
841 sage: n = ZZ.random_element(5)
842 sage: F = QuadraticField(-1, 'i')
843 sage: X = random_matrix(F, n)
844 sage: Y = random_matrix(F, n)
845 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
846 sage: expected = _embed_complex_matrix(X*Y)
847 sage: actual == expected
853 raise ValueError("the matrix 'M' must be square")
854 field
= M
.base_ring()
859 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
861 # We can drop the imaginaries here.
862 return matrix
.block(field
.base_ring(), n
, blocks
)
865 def _unembed_complex_matrix(M
):
867 The inverse of _embed_complex_matrix().
871 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
872 ....: _unembed_complex_matrix)
876 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
877 ....: [-2, 1, -4, 3],
878 ....: [ 9, 10, 11, 12],
879 ....: [-10, 9, -12, 11] ])
880 sage: _unembed_complex_matrix(A)
882 [ 10*i + 9 12*i + 11]
886 Unembedding is the inverse of embedding::
888 sage: set_random_seed()
889 sage: F = QuadraticField(-1, 'i')
890 sage: M = random_matrix(F, 3)
891 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
897 raise ValueError("the matrix 'M' must be square")
898 if not n
.mod(2).is_zero():
899 raise ValueError("the matrix 'M' must be a complex embedding")
901 F
= QuadraticField(-1, 'i')
904 # Go top-left to bottom-right (reading order), converting every
905 # 2-by-2 block we see to a single complex element.
907 for k
in xrange(n
/2):
908 for j
in xrange(n
/2):
909 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
910 if submat
[0,0] != submat
[1,1]:
911 raise ValueError('bad on-diagonal submatrix')
912 if submat
[0,1] != -submat
[1,0]:
913 raise ValueError('bad off-diagonal submatrix')
914 z
= submat
[0,0] + submat
[0,1]*i
917 return matrix(F
, n
/2, elements
)
920 def _embed_quaternion_matrix(M
):
922 Embed the n-by-n quaternion matrix ``M`` into the space of real
923 matrices of size 4n-by-4n by first sending each quaternion entry
924 `z = a + bi + cj + dk` to the block-complex matrix
925 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
930 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
934 sage: Q = QuaternionAlgebra(QQ,-1,-1)
935 sage: i,j,k = Q.gens()
936 sage: x = 1 + 2*i + 3*j + 4*k
937 sage: M = matrix(Q, 1, [[x]])
938 sage: _embed_quaternion_matrix(M)
944 Embedding is a homomorphism (isomorphism, in fact)::
946 sage: set_random_seed()
947 sage: n = ZZ.random_element(5)
948 sage: Q = QuaternionAlgebra(QQ,-1,-1)
949 sage: X = random_matrix(Q, n)
950 sage: Y = random_matrix(Q, n)
951 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
952 sage: expected = _embed_quaternion_matrix(X*Y)
953 sage: actual == expected
957 quaternions
= M
.base_ring()
960 raise ValueError("the matrix 'M' must be square")
962 F
= QuadraticField(-1, 'i')
967 t
= z
.coefficient_tuple()
972 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
973 [-c
+ d
*i
, a
- b
*i
]])
974 blocks
.append(_embed_complex_matrix(cplx_matrix
))
976 # We should have real entries by now, so use the realest field
977 # we've got for the return value.
978 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
981 def _unembed_quaternion_matrix(M
):
983 The inverse of _embed_quaternion_matrix().
987 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
988 ....: _unembed_quaternion_matrix)
992 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
993 ....: [-2, 1, -4, 3],
994 ....: [-3, 4, 1, -2],
995 ....: [-4, -3, 2, 1]])
996 sage: _unembed_quaternion_matrix(M)
997 [1 + 2*i + 3*j + 4*k]
1001 Unembedding is the inverse of embedding::
1003 sage: set_random_seed()
1004 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1005 sage: M = random_matrix(Q, 3)
1006 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1012 raise ValueError("the matrix 'M' must be square")
1013 if not n
.mod(4).is_zero():
1014 raise ValueError("the matrix 'M' must be a complex embedding")
1016 Q
= QuaternionAlgebra(QQ
,-1,-1)
1019 # Go top-left to bottom-right (reading order), converting every
1020 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1023 for l
in xrange(n
/4):
1024 for m
in xrange(n
/4):
1025 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1026 if submat
[0,0] != submat
[1,1].conjugate():
1027 raise ValueError('bad on-diagonal submatrix')
1028 if submat
[0,1] != -submat
[1,0].conjugate():
1029 raise ValueError('bad off-diagonal submatrix')
1030 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1031 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1034 return matrix(Q
, n
/4, elements
)
1037 # The usual inner product on R^n.
1039 return x
.to_vector().inner_product(y
.to_vector())
1041 # The inner product used for the real symmetric simple EJA.
1042 # We keep it as a separate function because e.g. the complex
1043 # algebra uses the same inner product, except divided by 2.
1044 def _matrix_ip(X
,Y
):
1045 X_mat
= X
.natural_representation()
1046 Y_mat
= Y
.natural_representation()
1047 return (X_mat
*Y_mat
).trace()
1050 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1052 The rank-n simple EJA consisting of real symmetric n-by-n
1053 matrices, the usual symmetric Jordan product, and the trace inner
1054 product. It has dimension `(n^2 + n)/2` over the reals.
1058 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1062 sage: J = RealSymmetricEJA(2)
1063 sage: e0, e1, e2 = J.gens()
1073 The dimension of this algebra is `(n^2 + n) / 2`::
1075 sage: set_random_seed()
1076 sage: n = ZZ.random_element(1,5)
1077 sage: J = RealSymmetricEJA(n)
1078 sage: J.dimension() == (n^2 + n)/2
1081 The Jordan multiplication is what we think it is::
1083 sage: set_random_seed()
1084 sage: n = ZZ.random_element(1,5)
1085 sage: J = RealSymmetricEJA(n)
1086 sage: x = J.random_element()
1087 sage: y = J.random_element()
1088 sage: actual = (x*y).natural_representation()
1089 sage: X = x.natural_representation()
1090 sage: Y = y.natural_representation()
1091 sage: expected = (X*Y + Y*X)/2
1092 sage: actual == expected
1094 sage: J(expected) == x*y
1098 def __init__(self
, n
, field
=QQ
):
1099 S
= _real_symmetric_basis(n
, field
=field
)
1100 Qs
= _multiplication_table_from_matrix_basis(S
)
1102 fdeja
= super(RealSymmetricEJA
, self
)
1103 return fdeja
.__init
__(field
,
1108 def inner_product(self
, x
, y
):
1109 return _matrix_ip(x
,y
)
1112 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1114 The rank-n simple EJA consisting of complex Hermitian n-by-n
1115 matrices over the real numbers, the usual symmetric Jordan product,
1116 and the real-part-of-trace inner product. It has dimension `n^2` over
1121 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1125 The dimension of this algebra is `n^2`::
1127 sage: set_random_seed()
1128 sage: n = ZZ.random_element(1,5)
1129 sage: J = ComplexHermitianEJA(n)
1130 sage: J.dimension() == n^2
1133 The Jordan multiplication is what we think it is::
1135 sage: set_random_seed()
1136 sage: n = ZZ.random_element(1,5)
1137 sage: J = ComplexHermitianEJA(n)
1138 sage: x = J.random_element()
1139 sage: y = J.random_element()
1140 sage: actual = (x*y).natural_representation()
1141 sage: X = x.natural_representation()
1142 sage: Y = y.natural_representation()
1143 sage: expected = (X*Y + Y*X)/2
1144 sage: actual == expected
1146 sage: J(expected) == x*y
1150 def __init__(self
, n
, field
=QQ
):
1151 S
= _complex_hermitian_basis(n
)
1152 Qs
= _multiplication_table_from_matrix_basis(S
)
1154 fdeja
= super(ComplexHermitianEJA
, self
)
1155 return fdeja
.__init
__(field
,
1161 def inner_product(self
, x
, y
):
1162 # Since a+bi on the diagonal is represented as
1167 # we'll double-count the "a" entries if we take the trace of
1169 return _matrix_ip(x
,y
)/2
1172 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1174 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1175 matrices, the usual symmetric Jordan product, and the
1176 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1181 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1185 The dimension of this algebra is `n^2`::
1187 sage: set_random_seed()
1188 sage: n = ZZ.random_element(1,5)
1189 sage: J = QuaternionHermitianEJA(n)
1190 sage: J.dimension() == 2*(n^2) - n
1193 The Jordan multiplication is what we think it is::
1195 sage: set_random_seed()
1196 sage: n = ZZ.random_element(1,5)
1197 sage: J = QuaternionHermitianEJA(n)
1198 sage: x = J.random_element()
1199 sage: y = J.random_element()
1200 sage: actual = (x*y).natural_representation()
1201 sage: X = x.natural_representation()
1202 sage: Y = y.natural_representation()
1203 sage: expected = (X*Y + Y*X)/2
1204 sage: actual == expected
1206 sage: J(expected) == x*y
1210 def __init__(self
, n
, field
=QQ
):
1211 S
= _quaternion_hermitian_basis(n
)
1212 Qs
= _multiplication_table_from_matrix_basis(S
)
1214 fdeja
= super(QuaternionHermitianEJA
, self
)
1215 return fdeja
.__init
__(field
,
1220 def inner_product(self
, x
, y
):
1221 # Since a+bi+cj+dk on the diagonal is represented as
1223 # a + bi +cj + dk = [ a b c d]
1228 # we'll quadruple-count the "a" entries if we take the trace of
1230 return _matrix_ip(x
,y
)/4
1233 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1235 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1236 with the usual inner product and jordan product ``x*y =
1237 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1242 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1246 This multiplication table can be verified by hand::
1248 sage: J = JordanSpinEJA(4)
1249 sage: e0,e1,e2,e3 = J.gens()
1266 def __init__(self
, n
, field
=QQ
):
1268 id_matrix
= matrix
.identity(field
, n
)
1270 ei
= id_matrix
.column(i
)
1271 Qi
= matrix
.zero(field
, n
)
1273 Qi
.set_column(0, ei
)
1274 Qi
+= matrix
.diagonal(n
, [ei
[0]]*n
)
1275 # The addition of the diagonal matrix adds an extra ei[0] in the
1276 # upper-left corner of the matrix.
1277 Qi
[0,0] = Qi
[0,0] * ~
field(2)
1280 # The rank of the spin algebra is two, unless we're in a
1281 # one-dimensional ambient space (because the rank is bounded by
1282 # the ambient dimension).
1283 fdeja
= super(JordanSpinEJA
, self
)
1284 return fdeja
.__init
__(field
, Qs
, rank
=min(n
,2))
1286 def inner_product(self
, x
, y
):
1287 return _usual_ip(x
,y
)