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
.magmatic_algebras
import MagmaticAlgebras
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
.misc
.table
import table
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
22 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
23 from mjo
.eja
.eja_utils
import _mat2vec
25 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
26 # This is an ugly hack needed to prevent the category framework
27 # from implementing a coercion from our base ring (e.g. the
28 # rationals) into the algebra. First of all -- such a coercion is
29 # nonsense to begin with. But more importantly, it tries to do so
30 # in the category of rings, and since our algebras aren't
31 # associative they generally won't be rings.
32 _no_generic_basering_coercion
= True
44 sage: from mjo.eja.eja_algebra import random_eja
48 By definition, Jordan multiplication commutes::
50 sage: set_random_seed()
51 sage: J = random_eja()
52 sage: x = J.random_element()
53 sage: y = J.random_element()
59 self
._natural
_basis
= natural_basis
62 category
= MagmaticAlgebras(field
).FiniteDimensional()
63 category
= category
.WithBasis().Unital()
65 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
67 range(len(mult_table
)),
70 self
.print_options(bracket
='')
72 # The multiplication table we're given is necessarily in terms
73 # of vectors, because we don't have an algebra yet for
74 # anything to be an element of. However, it's faster in the
75 # long run to have the multiplication table be in terms of
76 # algebra elements. We do this after calling the superclass
77 # constructor so that from_vector() knows what to do.
78 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
79 for ls
in mult_table
]
82 def _element_constructor_(self
, elt
):
84 Construct an element of this algebra from its natural
87 This gets called only after the parent element _call_ method
88 fails to find a coercion for the argument.
92 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
93 ....: RealCartesianProductEJA,
94 ....: RealSymmetricEJA)
98 The identity in `S^n` is converted to the identity in the EJA::
100 sage: J = RealSymmetricEJA(3)
101 sage: I = matrix.identity(QQ,3)
102 sage: J(I) == J.one()
105 This skew-symmetric matrix can't be represented in the EJA::
107 sage: J = RealSymmetricEJA(3)
108 sage: A = matrix(QQ,3, lambda i,j: i-j)
110 Traceback (most recent call last):
112 ArithmeticError: vector is not in free module
116 Ensure that we can convert any element of the two non-matrix
117 simple algebras (whose natural representations are their usual
118 vector representations) back and forth faithfully::
120 sage: set_random_seed()
121 sage: J = RealCartesianProductEJA(5)
122 sage: x = J.random_element()
123 sage: J(x.to_vector().column()) == x
125 sage: J = JordanSpinEJA(5)
126 sage: x = J.random_element()
127 sage: J(x.to_vector().column()) == x
132 # The superclass implementation of random_element()
133 # needs to be able to coerce "0" into the algebra.
136 natural_basis
= self
.natural_basis()
137 if elt
not in natural_basis
[0].matrix_space():
138 raise ValueError("not a naturally-represented algebra element")
140 # Thanks for nothing! Matrix spaces aren't vector
141 # spaces in Sage, so we have to figure out its
142 # natural-basis coordinates ourselves.
143 V
= VectorSpace(elt
.base_ring(), elt
.nrows()*elt
.ncols())
144 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
145 coords
= W
.coordinate_vector(_mat2vec(elt
))
146 return self
.from_vector(coords
)
151 Return a string representation of ``self``.
155 sage: from mjo.eja.eja_algebra import JordanSpinEJA
159 Ensure that it says what we think it says::
161 sage: JordanSpinEJA(2, field=QQ)
162 Euclidean Jordan algebra of degree 2 over Rational Field
163 sage: JordanSpinEJA(3, field=RDF)
164 Euclidean Jordan algebra of degree 3 over Real Double Field
167 # TODO: change this to say "dimension" and fix all the tests.
168 fmt
= "Euclidean Jordan algebra of degree {} over {}"
169 return fmt
.format(self
.dimension(), self
.base_ring())
171 def product_on_basis(self
, i
, j
):
172 return self
._multiplication
_table
[i
][j
]
174 def _a_regular_element(self
):
176 Guess a regular element. Needed to compute the basis for our
177 characteristic polynomial coefficients.
181 sage: from mjo.eja.eja_algebra import random_eja
185 Ensure that this hacky method succeeds for every algebra that we
186 know how to construct::
188 sage: set_random_seed()
189 sage: J = random_eja()
190 sage: J._a_regular_element().is_regular()
195 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
196 if not z
.is_regular():
197 raise ValueError("don't know a regular element")
202 def _charpoly_basis_space(self
):
204 Return the vector space spanned by the basis used in our
205 characteristic polynomial coefficients. This is used not only to
206 compute those coefficients, but also any time we need to
207 evaluate the coefficients (like when we compute the trace or
210 z
= self
._a
_regular
_element
()
211 V
= self
.vector_space()
212 V1
= V
.span_of_basis( (z
**k
).to_vector() for k
in range(self
.rank()) )
213 b
= (V1
.basis() + V1
.complement().basis())
214 return V
.span_of_basis(b
)
218 def _charpoly_coeff(self
, i
):
220 Return the coefficient polynomial "a_{i}" of this algebra's
221 general characteristic polynomial.
223 Having this be a separate cached method lets us compute and
224 store the trace/determinant (a_{r-1} and a_{0} respectively)
225 separate from the entire characteristic polynomial.
227 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
228 R
= A_of_x
.base_ring()
230 # Guaranteed by theory
233 # Danger: the in-place modification is done for performance
234 # reasons (reconstructing a matrix with huge polynomial
235 # entries is slow), but I don't know how cached_method works,
236 # so it's highly possible that we're modifying some global
237 # list variable by reference, here. In other words, you
238 # probably shouldn't call this method twice on the same
239 # algebra, at the same time, in two threads
240 Ai_orig
= A_of_x
.column(i
)
241 A_of_x
.set_column(i
,xr
)
242 numerator
= A_of_x
.det()
243 A_of_x
.set_column(i
,Ai_orig
)
245 # We're relying on the theory here to ensure that each a_i is
246 # indeed back in R, and the added negative signs are to make
247 # the whole charpoly expression sum to zero.
248 return R(-numerator
/detA
)
252 def _charpoly_matrix_system(self
):
254 Compute the matrix whose entries A_ij are polynomials in
255 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
256 corresponding to `x^r` and the determinent of the matrix A =
257 [A_ij]. In other words, all of the fixed (cachable) data needed
258 to compute the coefficients of the characteristic polynomial.
263 # Construct a new algebra over a multivariate polynomial ring...
264 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
265 R
= PolynomialRing(self
.base_ring(), names
)
266 # Hack around the fact that our multiplication table is in terms of
267 # algebra elements but the constructor wants it in terms of vectors.
268 vmt
= [ tuple(map(lambda x
: x
.to_vector(), ls
))
269 for ls
in self
._multiplication
_table
]
270 J
= FiniteDimensionalEuclideanJordanAlgebra(R
, tuple(vmt
), r
)
272 idmat
= matrix
.identity(J
.base_ring(), n
)
274 W
= self
._charpoly
_basis
_space
()
275 W
= W
.change_ring(R
.fraction_field())
277 # Starting with the standard coordinates x = (X1,X2,...,Xn)
278 # and then converting the entries to W-coordinates allows us
279 # to pass in the standard coordinates to the charpoly and get
280 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
283 # W.coordinates(x^2) eval'd at (standard z-coords)
287 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
289 # We want the middle equivalent thing in our matrix, but use
290 # the first equivalent thing instead so that we can pass in
291 # standard coordinates.
292 x
= J
.from_vector(W(R
.gens()))
294 # Handle the zeroth power separately, because computing
295 # the unit element in J is mathematically suspect.
296 x0
= W
.coordinate_vector(self
.one().to_vector())
298 l1
+= [ W
.coordinate_vector((x
**k
).to_vector()).column()
299 for k
in range(1,r
) ]
300 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
301 A_of_x
= matrix
.block(R
, 1, n
, (l1
+ l2
))
302 xr
= W
.coordinate_vector((x
**r
).to_vector())
303 return (A_of_x
, x
, xr
, A_of_x
.det())
307 def characteristic_polynomial(self
):
309 Return a characteristic polynomial that works for all elements
312 The resulting polynomial has `n+1` variables, where `n` is the
313 dimension of this algebra. The first `n` variables correspond to
314 the coordinates of an algebra element: when evaluated at the
315 coordinates of an algebra element with respect to a certain
316 basis, the result is a univariate polynomial (in the one
317 remaining variable ``t``), namely the characteristic polynomial
322 sage: from mjo.eja.eja_algebra import JordanSpinEJA
326 The characteristic polynomial in the spin algebra is given in
327 Alizadeh, Example 11.11::
329 sage: J = JordanSpinEJA(3)
330 sage: p = J.characteristic_polynomial(); p
331 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
332 sage: xvec = J.one().to_vector()
340 # The list of coefficient polynomials a_1, a_2, ..., a_n.
341 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
343 # We go to a bit of trouble here to reorder the
344 # indeterminates, so that it's easier to evaluate the
345 # characteristic polynomial at x's coordinates and get back
346 # something in terms of t, which is what we want.
348 S
= PolynomialRing(self
.base_ring(),'t')
350 S
= PolynomialRing(S
, R
.variable_names())
353 # Note: all entries past the rth should be zero. The
354 # coefficient of the highest power (x^r) is 1, but it doesn't
355 # appear in the solution vector which contains coefficients
356 # for the other powers (to make them sum to x^r).
358 a
[r
] = 1 # corresponds to x^r
360 # When the rank is equal to the dimension, trying to
361 # assign a[r] goes out-of-bounds.
362 a
.append(1) # corresponds to x^r
364 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
367 def inner_product(self
, x
, y
):
369 The inner product associated with this Euclidean Jordan algebra.
371 Defaults to the trace inner product, but can be overridden by
372 subclasses if they are sure that the necessary properties are
377 sage: from mjo.eja.eja_algebra import random_eja
381 The inner product must satisfy its axiom for this algebra to truly
382 be a Euclidean Jordan Algebra::
384 sage: set_random_seed()
385 sage: J = random_eja()
386 sage: x = J.random_element()
387 sage: y = J.random_element()
388 sage: z = J.random_element()
389 sage: (x*y).inner_product(z) == y.inner_product(x*z)
393 if (not x
in self
) or (not y
in self
):
394 raise TypeError("arguments must live in this algebra")
395 return x
.trace_inner_product(y
)
398 def multiplication_table(self
):
400 Return a visual representation of this algebra's multiplication
401 table (on basis elements).
405 sage: from mjo.eja.eja_algebra import JordanSpinEJA
409 sage: J = JordanSpinEJA(4)
410 sage: J.multiplication_table()
411 +----++----+----+----+----+
412 | * || e0 | e1 | e2 | e3 |
413 +====++====+====+====+====+
414 | e0 || e0 | e1 | e2 | e3 |
415 +----++----+----+----+----+
416 | e1 || e1 | e0 | 0 | 0 |
417 +----++----+----+----+----+
418 | e2 || e2 | 0 | e0 | 0 |
419 +----++----+----+----+----+
420 | e3 || e3 | 0 | 0 | e0 |
421 +----++----+----+----+----+
424 M
= list(self
._multiplication
_table
) # copy
425 for i
in range(len(M
)):
426 # M had better be "square"
427 M
[i
] = [self
.monomial(i
)] + M
[i
]
428 M
= [["*"] + list(self
.gens())] + M
429 return table(M
, header_row
=True, header_column
=True, frame
=True)
432 def natural_basis(self
):
434 Return a more-natural representation of this algebra's basis.
436 Every finite-dimensional Euclidean Jordan Algebra is a direct
437 sum of five simple algebras, four of which comprise Hermitian
438 matrices. This method returns the original "natural" basis
439 for our underlying vector space. (Typically, the natural basis
440 is used to construct the multiplication table in the first place.)
442 Note that this will always return a matrix. The standard basis
443 in `R^n` will be returned as `n`-by-`1` column matrices.
447 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
448 ....: RealSymmetricEJA)
452 sage: J = RealSymmetricEJA(2)
454 Finite family {0: e0, 1: e1, 2: e2}
455 sage: J.natural_basis()
463 sage: J = JordanSpinEJA(2)
465 Finite family {0: e0, 1: e1}
466 sage: J.natural_basis()
473 if self
._natural
_basis
is None:
474 return tuple( b
.to_vector().column() for b
in self
.basis() )
476 return self
._natural
_basis
482 Return the unit element of this algebra.
486 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
491 sage: J = RealCartesianProductEJA(5)
493 e0 + e1 + e2 + e3 + e4
497 The identity element acts like the identity::
499 sage: set_random_seed()
500 sage: J = random_eja()
501 sage: x = J.random_element()
502 sage: J.one()*x == x and x*J.one() == x
505 The matrix of the unit element's operator is the identity::
507 sage: set_random_seed()
508 sage: J = random_eja()
509 sage: actual = J.one().operator().matrix()
510 sage: expected = matrix.identity(J.base_ring(), J.dimension())
511 sage: actual == expected
515 # We can brute-force compute the matrices of the operators
516 # that correspond to the basis elements of this algebra.
517 # If some linear combination of those basis elements is the
518 # algebra identity, then the same linear combination of
519 # their matrices has to be the identity matrix.
521 # Of course, matrices aren't vectors in sage, so we have to
522 # appeal to the "long vectors" isometry.
523 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
525 # Now we use basis linear algebra to find the coefficients,
526 # of the matrices-as-vectors-linear-combination, which should
527 # work for the original algebra basis too.
528 A
= matrix
.column(self
.base_ring(), oper_vecs
)
530 # We used the isometry on the left-hand side already, but we
531 # still need to do it for the right-hand side. Recall that we
532 # wanted something that summed to the identity matrix.
533 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
535 # Now if there's an identity element in the algebra, this should work.
536 coeffs
= A
.solve_right(b
)
537 return self
.linear_combination(zip(self
.gens(), coeffs
))
542 Return the rank of this EJA.
546 The author knows of no algorithm to compute the rank of an EJA
547 where only the multiplication table is known. In lieu of one, we
548 require the rank to be specified when the algebra is created,
549 and simply pass along that number here.
553 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
554 ....: RealSymmetricEJA,
555 ....: ComplexHermitianEJA,
556 ....: QuaternionHermitianEJA,
561 The rank of the Jordan spin algebra is always two::
563 sage: JordanSpinEJA(2).rank()
565 sage: JordanSpinEJA(3).rank()
567 sage: JordanSpinEJA(4).rank()
570 The rank of the `n`-by-`n` Hermitian real, complex, or
571 quaternion matrices is `n`::
573 sage: RealSymmetricEJA(2).rank()
575 sage: ComplexHermitianEJA(2).rank()
577 sage: QuaternionHermitianEJA(2).rank()
579 sage: RealSymmetricEJA(5).rank()
581 sage: ComplexHermitianEJA(5).rank()
583 sage: QuaternionHermitianEJA(5).rank()
588 Ensure that every EJA that we know how to construct has a
589 positive integer rank::
591 sage: set_random_seed()
592 sage: r = random_eja().rank()
593 sage: r in ZZ and r > 0
600 def vector_space(self
):
602 Return the vector space that underlies this algebra.
606 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
610 sage: J = RealSymmetricEJA(2)
611 sage: J.vector_space()
612 Vector space of dimension 3 over Rational Field
615 return self
.zero().to_vector().parent().ambient_vector_space()
618 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
621 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
623 Return the Euclidean Jordan Algebra corresponding to the set
624 `R^n` under the Hadamard product.
626 Note: this is nothing more than the Cartesian product of ``n``
627 copies of the spin algebra. Once Cartesian product algebras
628 are implemented, this can go.
632 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
636 This multiplication table can be verified by hand::
638 sage: J = RealCartesianProductEJA(3)
639 sage: e0,e1,e2 = J.gens()
654 def __init__(self
, n
, field
=QQ
):
655 V
= VectorSpace(field
, n
)
656 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
659 fdeja
= super(RealCartesianProductEJA
, self
)
660 return fdeja
.__init
__(field
, mult_table
, rank
=n
)
662 def inner_product(self
, x
, y
):
663 return _usual_ip(x
,y
)
668 Return a "random" finite-dimensional Euclidean Jordan Algebra.
672 For now, we choose a random natural number ``n`` (greater than zero)
673 and then give you back one of the following:
675 * The cartesian product of the rational numbers ``n`` times; this is
676 ``QQ^n`` with the Hadamard product.
678 * The Jordan spin algebra on ``QQ^n``.
680 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
683 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
684 in the space of ``2n``-by-``2n`` real symmetric matrices.
686 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
687 in the space of ``4n``-by-``4n`` real symmetric matrices.
689 Later this might be extended to return Cartesian products of the
694 sage: from mjo.eja.eja_algebra import random_eja
699 Euclidean Jordan algebra of degree...
703 # The max_n component lets us choose different upper bounds on the
704 # value "n" that gets passed to the constructor. This is needed
705 # because e.g. R^{10} is reasonable to test, while the Hermitian
706 # 10-by-10 quaternion matrices are not.
707 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
709 (RealSymmetricEJA
, 5),
710 (ComplexHermitianEJA
, 4),
711 (QuaternionHermitianEJA
, 3)])
712 n
= ZZ
.random_element(1, max_n
)
713 return constructor(n
, field
=QQ
)
717 def _real_symmetric_basis(n
, field
=QQ
):
719 Return a basis for the space of real symmetric n-by-n matrices.
721 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
725 for j
in xrange(i
+1):
726 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
730 # Beware, orthogonal but not normalized!
731 Sij
= Eij
+ Eij
.transpose()
736 def _complex_hermitian_basis(n
, field
=QQ
):
738 Returns a basis for the space of complex Hermitian n-by-n matrices.
742 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
746 sage: set_random_seed()
747 sage: n = ZZ.random_element(1,5)
748 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
752 F
= QuadraticField(-1, 'I')
755 # This is like the symmetric case, but we need to be careful:
757 # * We want conjugate-symmetry, not just symmetry.
758 # * The diagonal will (as a result) be real.
762 for j
in xrange(i
+1):
763 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
765 Sij
= _embed_complex_matrix(Eij
)
768 # Beware, orthogonal but not normalized! The second one
769 # has a minus because it's conjugated.
770 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
772 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
777 def _quaternion_hermitian_basis(n
, field
=QQ
):
779 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
783 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
787 sage: set_random_seed()
788 sage: n = ZZ.random_element(1,5)
789 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
793 Q
= QuaternionAlgebra(QQ
,-1,-1)
796 # This is like the symmetric case, but we need to be careful:
798 # * We want conjugate-symmetry, not just symmetry.
799 # * The diagonal will (as a result) be real.
803 for j
in xrange(i
+1):
804 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
806 Sij
= _embed_quaternion_matrix(Eij
)
809 # Beware, orthogonal but not normalized! The second,
810 # third, and fourth ones have a minus because they're
812 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
814 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
816 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
818 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
824 def _multiplication_table_from_matrix_basis(basis
):
826 At least three of the five simple Euclidean Jordan algebras have the
827 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
828 multiplication on the right is matrix multiplication. Given a basis
829 for the underlying matrix space, this function returns a
830 multiplication table (obtained by looping through the basis
831 elements) for an algebra of those matrices.
833 # In S^2, for example, we nominally have four coordinates even
834 # though the space is of dimension three only. The vector space V
835 # is supposed to hold the entire long vector, and the subspace W
836 # of V will be spanned by the vectors that arise from symmetric
837 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
838 field
= basis
[0].base_ring()
839 dimension
= basis
[0].nrows()
841 V
= VectorSpace(field
, dimension
**2)
842 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
844 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
847 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
848 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
853 def _embed_complex_matrix(M
):
855 Embed the n-by-n complex matrix ``M`` into the space of real
856 matrices of size 2n-by-2n via the map the sends each entry `z = a +
857 bi` to the block matrix ``[[a,b],[-b,a]]``.
861 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
865 sage: F = QuadraticField(-1,'i')
866 sage: x1 = F(4 - 2*i)
867 sage: x2 = F(1 + 2*i)
870 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
871 sage: _embed_complex_matrix(M)
880 Embedding is a homomorphism (isomorphism, in fact)::
882 sage: set_random_seed()
883 sage: n = ZZ.random_element(5)
884 sage: F = QuadraticField(-1, 'i')
885 sage: X = random_matrix(F, n)
886 sage: Y = random_matrix(F, n)
887 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
888 sage: expected = _embed_complex_matrix(X*Y)
889 sage: actual == expected
895 raise ValueError("the matrix 'M' must be square")
896 field
= M
.base_ring()
901 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
903 # We can drop the imaginaries here.
904 return matrix
.block(field
.base_ring(), n
, blocks
)
907 def _unembed_complex_matrix(M
):
909 The inverse of _embed_complex_matrix().
913 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
914 ....: _unembed_complex_matrix)
918 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
919 ....: [-2, 1, -4, 3],
920 ....: [ 9, 10, 11, 12],
921 ....: [-10, 9, -12, 11] ])
922 sage: _unembed_complex_matrix(A)
924 [ 10*i + 9 12*i + 11]
928 Unembedding is the inverse of embedding::
930 sage: set_random_seed()
931 sage: F = QuadraticField(-1, 'i')
932 sage: M = random_matrix(F, 3)
933 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
939 raise ValueError("the matrix 'M' must be square")
940 if not n
.mod(2).is_zero():
941 raise ValueError("the matrix 'M' must be a complex embedding")
943 F
= QuadraticField(-1, 'i')
946 # Go top-left to bottom-right (reading order), converting every
947 # 2-by-2 block we see to a single complex element.
949 for k
in xrange(n
/2):
950 for j
in xrange(n
/2):
951 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
952 if submat
[0,0] != submat
[1,1]:
953 raise ValueError('bad on-diagonal submatrix')
954 if submat
[0,1] != -submat
[1,0]:
955 raise ValueError('bad off-diagonal submatrix')
956 z
= submat
[0,0] + submat
[0,1]*i
959 return matrix(F
, n
/2, elements
)
962 def _embed_quaternion_matrix(M
):
964 Embed the n-by-n quaternion matrix ``M`` into the space of real
965 matrices of size 4n-by-4n by first sending each quaternion entry
966 `z = a + bi + cj + dk` to the block-complex matrix
967 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
972 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
976 sage: Q = QuaternionAlgebra(QQ,-1,-1)
977 sage: i,j,k = Q.gens()
978 sage: x = 1 + 2*i + 3*j + 4*k
979 sage: M = matrix(Q, 1, [[x]])
980 sage: _embed_quaternion_matrix(M)
986 Embedding is a homomorphism (isomorphism, in fact)::
988 sage: set_random_seed()
989 sage: n = ZZ.random_element(5)
990 sage: Q = QuaternionAlgebra(QQ,-1,-1)
991 sage: X = random_matrix(Q, n)
992 sage: Y = random_matrix(Q, n)
993 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
994 sage: expected = _embed_quaternion_matrix(X*Y)
995 sage: actual == expected
999 quaternions
= M
.base_ring()
1002 raise ValueError("the matrix 'M' must be square")
1004 F
= QuadraticField(-1, 'i')
1009 t
= z
.coefficient_tuple()
1014 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1015 [-c
+ d
*i
, a
- b
*i
]])
1016 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1018 # We should have real entries by now, so use the realest field
1019 # we've got for the return value.
1020 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1023 def _unembed_quaternion_matrix(M
):
1025 The inverse of _embed_quaternion_matrix().
1029 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1030 ....: _unembed_quaternion_matrix)
1034 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1035 ....: [-2, 1, -4, 3],
1036 ....: [-3, 4, 1, -2],
1037 ....: [-4, -3, 2, 1]])
1038 sage: _unembed_quaternion_matrix(M)
1039 [1 + 2*i + 3*j + 4*k]
1043 Unembedding is the inverse of embedding::
1045 sage: set_random_seed()
1046 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1047 sage: M = random_matrix(Q, 3)
1048 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1054 raise ValueError("the matrix 'M' must be square")
1055 if not n
.mod(4).is_zero():
1056 raise ValueError("the matrix 'M' must be a complex embedding")
1058 Q
= QuaternionAlgebra(QQ
,-1,-1)
1061 # Go top-left to bottom-right (reading order), converting every
1062 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1065 for l
in xrange(n
/4):
1066 for m
in xrange(n
/4):
1067 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1068 if submat
[0,0] != submat
[1,1].conjugate():
1069 raise ValueError('bad on-diagonal submatrix')
1070 if submat
[0,1] != -submat
[1,0].conjugate():
1071 raise ValueError('bad off-diagonal submatrix')
1072 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1073 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1076 return matrix(Q
, n
/4, elements
)
1079 # The usual inner product on R^n.
1081 return x
.to_vector().inner_product(y
.to_vector())
1083 # The inner product used for the real symmetric simple EJA.
1084 # We keep it as a separate function because e.g. the complex
1085 # algebra uses the same inner product, except divided by 2.
1086 def _matrix_ip(X
,Y
):
1087 X_mat
= X
.natural_representation()
1088 Y_mat
= Y
.natural_representation()
1089 return (X_mat
*Y_mat
).trace()
1092 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1094 The rank-n simple EJA consisting of real symmetric n-by-n
1095 matrices, the usual symmetric Jordan product, and the trace inner
1096 product. It has dimension `(n^2 + n)/2` over the reals.
1100 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1104 sage: J = RealSymmetricEJA(2)
1105 sage: e0, e1, e2 = J.gens()
1115 The dimension of this algebra is `(n^2 + n) / 2`::
1117 sage: set_random_seed()
1118 sage: n = ZZ.random_element(1,5)
1119 sage: J = RealSymmetricEJA(n)
1120 sage: J.dimension() == (n^2 + n)/2
1123 The Jordan multiplication is what we think it is::
1125 sage: set_random_seed()
1126 sage: n = ZZ.random_element(1,5)
1127 sage: J = RealSymmetricEJA(n)
1128 sage: x = J.random_element()
1129 sage: y = J.random_element()
1130 sage: actual = (x*y).natural_representation()
1131 sage: X = x.natural_representation()
1132 sage: Y = y.natural_representation()
1133 sage: expected = (X*Y + Y*X)/2
1134 sage: actual == expected
1136 sage: J(expected) == x*y
1140 def __init__(self
, n
, field
=QQ
):
1141 S
= _real_symmetric_basis(n
, field
=field
)
1142 Qs
= _multiplication_table_from_matrix_basis(S
)
1144 fdeja
= super(RealSymmetricEJA
, self
)
1145 return fdeja
.__init
__(field
,
1150 def inner_product(self
, x
, y
):
1151 return _matrix_ip(x
,y
)
1154 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1156 The rank-n simple EJA consisting of complex Hermitian n-by-n
1157 matrices over the real numbers, the usual symmetric Jordan product,
1158 and the real-part-of-trace inner product. It has dimension `n^2` over
1163 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1167 The dimension of this algebra is `n^2`::
1169 sage: set_random_seed()
1170 sage: n = ZZ.random_element(1,5)
1171 sage: J = ComplexHermitianEJA(n)
1172 sage: J.dimension() == n^2
1175 The Jordan multiplication is what we think it is::
1177 sage: set_random_seed()
1178 sage: n = ZZ.random_element(1,5)
1179 sage: J = ComplexHermitianEJA(n)
1180 sage: x = J.random_element()
1181 sage: y = J.random_element()
1182 sage: actual = (x*y).natural_representation()
1183 sage: X = x.natural_representation()
1184 sage: Y = y.natural_representation()
1185 sage: expected = (X*Y + Y*X)/2
1186 sage: actual == expected
1188 sage: J(expected) == x*y
1192 def __init__(self
, n
, field
=QQ
):
1193 S
= _complex_hermitian_basis(n
)
1194 Qs
= _multiplication_table_from_matrix_basis(S
)
1196 fdeja
= super(ComplexHermitianEJA
, self
)
1197 return fdeja
.__init
__(field
,
1203 def inner_product(self
, x
, y
):
1204 # Since a+bi on the diagonal is represented as
1209 # we'll double-count the "a" entries if we take the trace of
1211 return _matrix_ip(x
,y
)/2
1214 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1216 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1217 matrices, the usual symmetric Jordan product, and the
1218 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1223 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1227 The dimension of this algebra is `n^2`::
1229 sage: set_random_seed()
1230 sage: n = ZZ.random_element(1,5)
1231 sage: J = QuaternionHermitianEJA(n)
1232 sage: J.dimension() == 2*(n^2) - n
1235 The Jordan multiplication is what we think it is::
1237 sage: set_random_seed()
1238 sage: n = ZZ.random_element(1,5)
1239 sage: J = QuaternionHermitianEJA(n)
1240 sage: x = J.random_element()
1241 sage: y = J.random_element()
1242 sage: actual = (x*y).natural_representation()
1243 sage: X = x.natural_representation()
1244 sage: Y = y.natural_representation()
1245 sage: expected = (X*Y + Y*X)/2
1246 sage: actual == expected
1248 sage: J(expected) == x*y
1252 def __init__(self
, n
, field
=QQ
):
1253 S
= _quaternion_hermitian_basis(n
)
1254 Qs
= _multiplication_table_from_matrix_basis(S
)
1256 fdeja
= super(QuaternionHermitianEJA
, self
)
1257 return fdeja
.__init
__(field
,
1262 def inner_product(self
, x
, y
):
1263 # Since a+bi+cj+dk on the diagonal is represented as
1265 # a + bi +cj + dk = [ a b c d]
1270 # we'll quadruple-count the "a" entries if we take the trace of
1272 return _matrix_ip(x
,y
)/4
1275 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1277 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1278 with the usual inner product and jordan product ``x*y =
1279 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1284 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1288 This multiplication table can be verified by hand::
1290 sage: J = JordanSpinEJA(4)
1291 sage: e0,e1,e2,e3 = J.gens()
1308 def __init__(self
, n
, field
=QQ
):
1309 V
= VectorSpace(field
, n
)
1310 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1320 z0
= x
.inner_product(y
)
1321 zbar
= y0
*xbar
+ x0
*ybar
1322 z
= V([z0
] + zbar
.list())
1323 mult_table
[i
][j
] = z
1325 # The rank of the spin algebra is two, unless we're in a
1326 # one-dimensional ambient space (because the rank is bounded by
1327 # the ambient dimension).
1328 fdeja
= super(JordanSpinEJA
, self
)
1329 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2))
1331 def inner_product(self
, x
, y
):
1332 return _usual_ip(x
,y
)