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 dimension 2 over Rational Field
163 sage: JordanSpinEJA(3, field=RDF)
164 Euclidean Jordan algebra of dimension 3 over Real Double Field
167 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
168 return fmt
.format(self
.dimension(), self
.base_ring())
170 def product_on_basis(self
, i
, j
):
171 return self
._multiplication
_table
[i
][j
]
173 def _a_regular_element(self
):
175 Guess a regular element. Needed to compute the basis for our
176 characteristic polynomial coefficients.
180 sage: from mjo.eja.eja_algebra import random_eja
184 Ensure that this hacky method succeeds for every algebra that we
185 know how to construct::
187 sage: set_random_seed()
188 sage: J = random_eja()
189 sage: J._a_regular_element().is_regular()
194 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
195 if not z
.is_regular():
196 raise ValueError("don't know a regular element")
201 def _charpoly_basis_space(self
):
203 Return the vector space spanned by the basis used in our
204 characteristic polynomial coefficients. This is used not only to
205 compute those coefficients, but also any time we need to
206 evaluate the coefficients (like when we compute the trace or
209 z
= self
._a
_regular
_element
()
210 V
= self
.vector_space()
211 V1
= V
.span_of_basis( (z
**k
).to_vector() for k
in range(self
.rank()) )
212 b
= (V1
.basis() + V1
.complement().basis())
213 return V
.span_of_basis(b
)
217 def _charpoly_coeff(self
, i
):
219 Return the coefficient polynomial "a_{i}" of this algebra's
220 general characteristic polynomial.
222 Having this be a separate cached method lets us compute and
223 store the trace/determinant (a_{r-1} and a_{0} respectively)
224 separate from the entire characteristic polynomial.
226 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
227 R
= A_of_x
.base_ring()
229 # Guaranteed by theory
232 # Danger: the in-place modification is done for performance
233 # reasons (reconstructing a matrix with huge polynomial
234 # entries is slow), but I don't know how cached_method works,
235 # so it's highly possible that we're modifying some global
236 # list variable by reference, here. In other words, you
237 # probably shouldn't call this method twice on the same
238 # algebra, at the same time, in two threads
239 Ai_orig
= A_of_x
.column(i
)
240 A_of_x
.set_column(i
,xr
)
241 numerator
= A_of_x
.det()
242 A_of_x
.set_column(i
,Ai_orig
)
244 # We're relying on the theory here to ensure that each a_i is
245 # indeed back in R, and the added negative signs are to make
246 # the whole charpoly expression sum to zero.
247 return R(-numerator
/detA
)
251 def _charpoly_matrix_system(self
):
253 Compute the matrix whose entries A_ij are polynomials in
254 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
255 corresponding to `x^r` and the determinent of the matrix A =
256 [A_ij]. In other words, all of the fixed (cachable) data needed
257 to compute the coefficients of the characteristic polynomial.
262 # Construct a new algebra over a multivariate polynomial ring...
263 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
264 R
= PolynomialRing(self
.base_ring(), names
)
265 # Hack around the fact that our multiplication table is in terms of
266 # algebra elements but the constructor wants it in terms of vectors.
267 vmt
= [ tuple(map(lambda x
: x
.to_vector(), ls
))
268 for ls
in self
._multiplication
_table
]
269 J
= FiniteDimensionalEuclideanJordanAlgebra(R
, tuple(vmt
), r
)
271 idmat
= matrix
.identity(J
.base_ring(), n
)
273 W
= self
._charpoly
_basis
_space
()
274 W
= W
.change_ring(R
.fraction_field())
276 # Starting with the standard coordinates x = (X1,X2,...,Xn)
277 # and then converting the entries to W-coordinates allows us
278 # to pass in the standard coordinates to the charpoly and get
279 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
282 # W.coordinates(x^2) eval'd at (standard z-coords)
286 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
288 # We want the middle equivalent thing in our matrix, but use
289 # the first equivalent thing instead so that we can pass in
290 # standard coordinates.
291 x
= J
.from_vector(W(R
.gens()))
293 # Handle the zeroth power separately, because computing
294 # the unit element in J is mathematically suspect.
295 x0
= W
.coordinate_vector(self
.one().to_vector())
297 l1
+= [ W
.coordinate_vector((x
**k
).to_vector()).column()
298 for k
in range(1,r
) ]
299 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
300 A_of_x
= matrix
.block(R
, 1, n
, (l1
+ l2
))
301 xr
= W
.coordinate_vector((x
**r
).to_vector())
302 return (A_of_x
, x
, xr
, A_of_x
.det())
306 def characteristic_polynomial(self
):
308 Return a characteristic polynomial that works for all elements
311 The resulting polynomial has `n+1` variables, where `n` is the
312 dimension of this algebra. The first `n` variables correspond to
313 the coordinates of an algebra element: when evaluated at the
314 coordinates of an algebra element with respect to a certain
315 basis, the result is a univariate polynomial (in the one
316 remaining variable ``t``), namely the characteristic polynomial
321 sage: from mjo.eja.eja_algebra import JordanSpinEJA
325 The characteristic polynomial in the spin algebra is given in
326 Alizadeh, Example 11.11::
328 sage: J = JordanSpinEJA(3)
329 sage: p = J.characteristic_polynomial(); p
330 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
331 sage: xvec = J.one().to_vector()
339 # The list of coefficient polynomials a_1, a_2, ..., a_n.
340 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
342 # We go to a bit of trouble here to reorder the
343 # indeterminates, so that it's easier to evaluate the
344 # characteristic polynomial at x's coordinates and get back
345 # something in terms of t, which is what we want.
347 S
= PolynomialRing(self
.base_ring(),'t')
349 S
= PolynomialRing(S
, R
.variable_names())
352 # Note: all entries past the rth should be zero. The
353 # coefficient of the highest power (x^r) is 1, but it doesn't
354 # appear in the solution vector which contains coefficients
355 # for the other powers (to make them sum to x^r).
357 a
[r
] = 1 # corresponds to x^r
359 # When the rank is equal to the dimension, trying to
360 # assign a[r] goes out-of-bounds.
361 a
.append(1) # corresponds to x^r
363 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
366 def inner_product(self
, x
, y
):
368 The inner product associated with this Euclidean Jordan algebra.
370 Defaults to the trace inner product, but can be overridden by
371 subclasses if they are sure that the necessary properties are
376 sage: from mjo.eja.eja_algebra import random_eja
380 The inner product must satisfy its axiom for this algebra to truly
381 be a Euclidean Jordan Algebra::
383 sage: set_random_seed()
384 sage: J = random_eja()
385 sage: x = J.random_element()
386 sage: y = J.random_element()
387 sage: z = J.random_element()
388 sage: (x*y).inner_product(z) == y.inner_product(x*z)
392 if (not x
in self
) or (not y
in self
):
393 raise TypeError("arguments must live in this algebra")
394 return x
.trace_inner_product(y
)
397 def multiplication_table(self
):
399 Return a visual representation of this algebra's multiplication
400 table (on basis elements).
404 sage: from mjo.eja.eja_algebra import JordanSpinEJA
408 sage: J = JordanSpinEJA(4)
409 sage: J.multiplication_table()
410 +----++----+----+----+----+
411 | * || e0 | e1 | e2 | e3 |
412 +====++====+====+====+====+
413 | e0 || e0 | e1 | e2 | e3 |
414 +----++----+----+----+----+
415 | e1 || e1 | e0 | 0 | 0 |
416 +----++----+----+----+----+
417 | e2 || e2 | 0 | e0 | 0 |
418 +----++----+----+----+----+
419 | e3 || e3 | 0 | 0 | e0 |
420 +----++----+----+----+----+
423 M
= list(self
._multiplication
_table
) # copy
424 for i
in range(len(M
)):
425 # M had better be "square"
426 M
[i
] = [self
.monomial(i
)] + M
[i
]
427 M
= [["*"] + list(self
.gens())] + M
428 return table(M
, header_row
=True, header_column
=True, frame
=True)
431 def natural_basis(self
):
433 Return a more-natural representation of this algebra's basis.
435 Every finite-dimensional Euclidean Jordan Algebra is a direct
436 sum of five simple algebras, four of which comprise Hermitian
437 matrices. This method returns the original "natural" basis
438 for our underlying vector space. (Typically, the natural basis
439 is used to construct the multiplication table in the first place.)
441 Note that this will always return a matrix. The standard basis
442 in `R^n` will be returned as `n`-by-`1` column matrices.
446 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
447 ....: RealSymmetricEJA)
451 sage: J = RealSymmetricEJA(2)
453 Finite family {0: e0, 1: e1, 2: e2}
454 sage: J.natural_basis()
462 sage: J = JordanSpinEJA(2)
464 Finite family {0: e0, 1: e1}
465 sage: J.natural_basis()
472 if self
._natural
_basis
is None:
473 return tuple( b
.to_vector().column() for b
in self
.basis() )
475 return self
._natural
_basis
481 Return the unit element of this algebra.
485 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
490 sage: J = RealCartesianProductEJA(5)
492 e0 + e1 + e2 + e3 + e4
496 The identity element acts like the identity::
498 sage: set_random_seed()
499 sage: J = random_eja()
500 sage: x = J.random_element()
501 sage: J.one()*x == x and x*J.one() == x
504 The matrix of the unit element's operator is the identity::
506 sage: set_random_seed()
507 sage: J = random_eja()
508 sage: actual = J.one().operator().matrix()
509 sage: expected = matrix.identity(J.base_ring(), J.dimension())
510 sage: actual == expected
514 # We can brute-force compute the matrices of the operators
515 # that correspond to the basis elements of this algebra.
516 # If some linear combination of those basis elements is the
517 # algebra identity, then the same linear combination of
518 # their matrices has to be the identity matrix.
520 # Of course, matrices aren't vectors in sage, so we have to
521 # appeal to the "long vectors" isometry.
522 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
524 # Now we use basis linear algebra to find the coefficients,
525 # of the matrices-as-vectors-linear-combination, which should
526 # work for the original algebra basis too.
527 A
= matrix
.column(self
.base_ring(), oper_vecs
)
529 # We used the isometry on the left-hand side already, but we
530 # still need to do it for the right-hand side. Recall that we
531 # wanted something that summed to the identity matrix.
532 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
534 # Now if there's an identity element in the algebra, this should work.
535 coeffs
= A
.solve_right(b
)
536 return self
.linear_combination(zip(self
.gens(), coeffs
))
541 Return the rank of this EJA.
545 The author knows of no algorithm to compute the rank of an EJA
546 where only the multiplication table is known. In lieu of one, we
547 require the rank to be specified when the algebra is created,
548 and simply pass along that number here.
552 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
553 ....: RealSymmetricEJA,
554 ....: ComplexHermitianEJA,
555 ....: QuaternionHermitianEJA,
560 The rank of the Jordan spin algebra is always two::
562 sage: JordanSpinEJA(2).rank()
564 sage: JordanSpinEJA(3).rank()
566 sage: JordanSpinEJA(4).rank()
569 The rank of the `n`-by-`n` Hermitian real, complex, or
570 quaternion matrices is `n`::
572 sage: RealSymmetricEJA(2).rank()
574 sage: ComplexHermitianEJA(2).rank()
576 sage: QuaternionHermitianEJA(2).rank()
578 sage: RealSymmetricEJA(5).rank()
580 sage: ComplexHermitianEJA(5).rank()
582 sage: QuaternionHermitianEJA(5).rank()
587 Ensure that every EJA that we know how to construct has a
588 positive integer rank::
590 sage: set_random_seed()
591 sage: r = random_eja().rank()
592 sage: r in ZZ and r > 0
599 def vector_space(self
):
601 Return the vector space that underlies this algebra.
605 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
609 sage: J = RealSymmetricEJA(2)
610 sage: J.vector_space()
611 Vector space of dimension 3 over Rational Field
614 return self
.zero().to_vector().parent().ambient_vector_space()
617 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
620 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
622 Return the Euclidean Jordan Algebra corresponding to the set
623 `R^n` under the Hadamard product.
625 Note: this is nothing more than the Cartesian product of ``n``
626 copies of the spin algebra. Once Cartesian product algebras
627 are implemented, this can go.
631 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
635 This multiplication table can be verified by hand::
637 sage: J = RealCartesianProductEJA(3)
638 sage: e0,e1,e2 = J.gens()
653 def __init__(self
, n
, field
=QQ
):
654 V
= VectorSpace(field
, n
)
655 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
658 fdeja
= super(RealCartesianProductEJA
, self
)
659 return fdeja
.__init
__(field
, mult_table
, rank
=n
)
661 def inner_product(self
, x
, y
):
662 return _usual_ip(x
,y
)
667 Return a "random" finite-dimensional Euclidean Jordan Algebra.
671 For now, we choose a random natural number ``n`` (greater than zero)
672 and then give you back one of the following:
674 * The cartesian product of the rational numbers ``n`` times; this is
675 ``QQ^n`` with the Hadamard product.
677 * The Jordan spin algebra on ``QQ^n``.
679 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
682 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
683 in the space of ``2n``-by-``2n`` real symmetric matrices.
685 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
686 in the space of ``4n``-by-``4n`` real symmetric matrices.
688 Later this might be extended to return Cartesian products of the
693 sage: from mjo.eja.eja_algebra import random_eja
698 Euclidean Jordan algebra of dimension...
702 # The max_n component lets us choose different upper bounds on the
703 # value "n" that gets passed to the constructor. This is needed
704 # because e.g. R^{10} is reasonable to test, while the Hermitian
705 # 10-by-10 quaternion matrices are not.
706 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
708 (RealSymmetricEJA
, 5),
709 (ComplexHermitianEJA
, 4),
710 (QuaternionHermitianEJA
, 3)])
711 n
= ZZ
.random_element(1, max_n
)
712 return constructor(n
, field
=QQ
)
716 def _real_symmetric_basis(n
, field
=QQ
):
718 Return a basis for the space of real symmetric n-by-n matrices.
720 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
724 for j
in xrange(i
+1):
725 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
729 # Beware, orthogonal but not normalized!
730 Sij
= Eij
+ Eij
.transpose()
735 def _complex_hermitian_basis(n
, field
=QQ
):
737 Returns a basis for the space of complex Hermitian n-by-n matrices.
741 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
745 sage: set_random_seed()
746 sage: n = ZZ.random_element(1,5)
747 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
751 F
= QuadraticField(-1, 'I')
754 # This is like the symmetric case, but we need to be careful:
756 # * We want conjugate-symmetry, not just symmetry.
757 # * The diagonal will (as a result) be real.
761 for j
in xrange(i
+1):
762 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
764 Sij
= _embed_complex_matrix(Eij
)
767 # Beware, orthogonal but not normalized! The second one
768 # has a minus because it's conjugated.
769 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
771 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
776 def _quaternion_hermitian_basis(n
, field
=QQ
):
778 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
782 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
786 sage: set_random_seed()
787 sage: n = ZZ.random_element(1,5)
788 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
792 Q
= QuaternionAlgebra(QQ
,-1,-1)
795 # This is like the symmetric case, but we need to be careful:
797 # * We want conjugate-symmetry, not just symmetry.
798 # * The diagonal will (as a result) be real.
802 for j
in xrange(i
+1):
803 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
805 Sij
= _embed_quaternion_matrix(Eij
)
808 # Beware, orthogonal but not normalized! The second,
809 # third, and fourth ones have a minus because they're
811 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
813 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
815 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
817 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
823 def _multiplication_table_from_matrix_basis(basis
):
825 At least three of the five simple Euclidean Jordan algebras have the
826 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
827 multiplication on the right is matrix multiplication. Given a basis
828 for the underlying matrix space, this function returns a
829 multiplication table (obtained by looping through the basis
830 elements) for an algebra of those matrices.
832 # In S^2, for example, we nominally have four coordinates even
833 # though the space is of dimension three only. The vector space V
834 # is supposed to hold the entire long vector, and the subspace W
835 # of V will be spanned by the vectors that arise from symmetric
836 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
837 field
= basis
[0].base_ring()
838 dimension
= basis
[0].nrows()
840 V
= VectorSpace(field
, dimension
**2)
841 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
843 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
846 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
847 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
852 def _embed_complex_matrix(M
):
854 Embed the n-by-n complex matrix ``M`` into the space of real
855 matrices of size 2n-by-2n via the map the sends each entry `z = a +
856 bi` to the block matrix ``[[a,b],[-b,a]]``.
860 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
864 sage: F = QuadraticField(-1,'i')
865 sage: x1 = F(4 - 2*i)
866 sage: x2 = F(1 + 2*i)
869 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
870 sage: _embed_complex_matrix(M)
879 Embedding is a homomorphism (isomorphism, in fact)::
881 sage: set_random_seed()
882 sage: n = ZZ.random_element(5)
883 sage: F = QuadraticField(-1, 'i')
884 sage: X = random_matrix(F, n)
885 sage: Y = random_matrix(F, n)
886 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
887 sage: expected = _embed_complex_matrix(X*Y)
888 sage: actual == expected
894 raise ValueError("the matrix 'M' must be square")
895 field
= M
.base_ring()
900 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
902 # We can drop the imaginaries here.
903 return matrix
.block(field
.base_ring(), n
, blocks
)
906 def _unembed_complex_matrix(M
):
908 The inverse of _embed_complex_matrix().
912 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
913 ....: _unembed_complex_matrix)
917 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
918 ....: [-2, 1, -4, 3],
919 ....: [ 9, 10, 11, 12],
920 ....: [-10, 9, -12, 11] ])
921 sage: _unembed_complex_matrix(A)
923 [ 10*i + 9 12*i + 11]
927 Unembedding is the inverse of embedding::
929 sage: set_random_seed()
930 sage: F = QuadraticField(-1, 'i')
931 sage: M = random_matrix(F, 3)
932 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
938 raise ValueError("the matrix 'M' must be square")
939 if not n
.mod(2).is_zero():
940 raise ValueError("the matrix 'M' must be a complex embedding")
942 F
= QuadraticField(-1, 'i')
945 # Go top-left to bottom-right (reading order), converting every
946 # 2-by-2 block we see to a single complex element.
948 for k
in xrange(n
/2):
949 for j
in xrange(n
/2):
950 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
951 if submat
[0,0] != submat
[1,1]:
952 raise ValueError('bad on-diagonal submatrix')
953 if submat
[0,1] != -submat
[1,0]:
954 raise ValueError('bad off-diagonal submatrix')
955 z
= submat
[0,0] + submat
[0,1]*i
958 return matrix(F
, n
/2, elements
)
961 def _embed_quaternion_matrix(M
):
963 Embed the n-by-n quaternion matrix ``M`` into the space of real
964 matrices of size 4n-by-4n by first sending each quaternion entry
965 `z = a + bi + cj + dk` to the block-complex matrix
966 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
971 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
975 sage: Q = QuaternionAlgebra(QQ,-1,-1)
976 sage: i,j,k = Q.gens()
977 sage: x = 1 + 2*i + 3*j + 4*k
978 sage: M = matrix(Q, 1, [[x]])
979 sage: _embed_quaternion_matrix(M)
985 Embedding is a homomorphism (isomorphism, in fact)::
987 sage: set_random_seed()
988 sage: n = ZZ.random_element(5)
989 sage: Q = QuaternionAlgebra(QQ,-1,-1)
990 sage: X = random_matrix(Q, n)
991 sage: Y = random_matrix(Q, n)
992 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
993 sage: expected = _embed_quaternion_matrix(X*Y)
994 sage: actual == expected
998 quaternions
= M
.base_ring()
1001 raise ValueError("the matrix 'M' must be square")
1003 F
= QuadraticField(-1, 'i')
1008 t
= z
.coefficient_tuple()
1013 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1014 [-c
+ d
*i
, a
- b
*i
]])
1015 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1017 # We should have real entries by now, so use the realest field
1018 # we've got for the return value.
1019 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1022 def _unembed_quaternion_matrix(M
):
1024 The inverse of _embed_quaternion_matrix().
1028 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1029 ....: _unembed_quaternion_matrix)
1033 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1034 ....: [-2, 1, -4, 3],
1035 ....: [-3, 4, 1, -2],
1036 ....: [-4, -3, 2, 1]])
1037 sage: _unembed_quaternion_matrix(M)
1038 [1 + 2*i + 3*j + 4*k]
1042 Unembedding is the inverse of embedding::
1044 sage: set_random_seed()
1045 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1046 sage: M = random_matrix(Q, 3)
1047 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1053 raise ValueError("the matrix 'M' must be square")
1054 if not n
.mod(4).is_zero():
1055 raise ValueError("the matrix 'M' must be a complex embedding")
1057 Q
= QuaternionAlgebra(QQ
,-1,-1)
1060 # Go top-left to bottom-right (reading order), converting every
1061 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1064 for l
in xrange(n
/4):
1065 for m
in xrange(n
/4):
1066 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1067 if submat
[0,0] != submat
[1,1].conjugate():
1068 raise ValueError('bad on-diagonal submatrix')
1069 if submat
[0,1] != -submat
[1,0].conjugate():
1070 raise ValueError('bad off-diagonal submatrix')
1071 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1072 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1075 return matrix(Q
, n
/4, elements
)
1078 # The usual inner product on R^n.
1080 return x
.to_vector().inner_product(y
.to_vector())
1082 # The inner product used for the real symmetric simple EJA.
1083 # We keep it as a separate function because e.g. the complex
1084 # algebra uses the same inner product, except divided by 2.
1085 def _matrix_ip(X
,Y
):
1086 X_mat
= X
.natural_representation()
1087 Y_mat
= Y
.natural_representation()
1088 return (X_mat
*Y_mat
).trace()
1091 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1093 The rank-n simple EJA consisting of real symmetric n-by-n
1094 matrices, the usual symmetric Jordan product, and the trace inner
1095 product. It has dimension `(n^2 + n)/2` over the reals.
1099 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1103 sage: J = RealSymmetricEJA(2)
1104 sage: e0, e1, e2 = J.gens()
1114 The dimension of this algebra is `(n^2 + n) / 2`::
1116 sage: set_random_seed()
1117 sage: n = ZZ.random_element(1,5)
1118 sage: J = RealSymmetricEJA(n)
1119 sage: J.dimension() == (n^2 + n)/2
1122 The Jordan multiplication is what we think it is::
1124 sage: set_random_seed()
1125 sage: n = ZZ.random_element(1,5)
1126 sage: J = RealSymmetricEJA(n)
1127 sage: x = J.random_element()
1128 sage: y = J.random_element()
1129 sage: actual = (x*y).natural_representation()
1130 sage: X = x.natural_representation()
1131 sage: Y = y.natural_representation()
1132 sage: expected = (X*Y + Y*X)/2
1133 sage: actual == expected
1135 sage: J(expected) == x*y
1139 def __init__(self
, n
, field
=QQ
):
1140 S
= _real_symmetric_basis(n
, field
=field
)
1141 Qs
= _multiplication_table_from_matrix_basis(S
)
1143 fdeja
= super(RealSymmetricEJA
, self
)
1144 return fdeja
.__init
__(field
,
1149 def inner_product(self
, x
, y
):
1150 return _matrix_ip(x
,y
)
1153 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1155 The rank-n simple EJA consisting of complex Hermitian n-by-n
1156 matrices over the real numbers, the usual symmetric Jordan product,
1157 and the real-part-of-trace inner product. It has dimension `n^2` over
1162 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1166 The dimension of this algebra is `n^2`::
1168 sage: set_random_seed()
1169 sage: n = ZZ.random_element(1,5)
1170 sage: J = ComplexHermitianEJA(n)
1171 sage: J.dimension() == n^2
1174 The Jordan multiplication is what we think it is::
1176 sage: set_random_seed()
1177 sage: n = ZZ.random_element(1,5)
1178 sage: J = ComplexHermitianEJA(n)
1179 sage: x = J.random_element()
1180 sage: y = J.random_element()
1181 sage: actual = (x*y).natural_representation()
1182 sage: X = x.natural_representation()
1183 sage: Y = y.natural_representation()
1184 sage: expected = (X*Y + Y*X)/2
1185 sage: actual == expected
1187 sage: J(expected) == x*y
1191 def __init__(self
, n
, field
=QQ
):
1192 S
= _complex_hermitian_basis(n
)
1193 Qs
= _multiplication_table_from_matrix_basis(S
)
1195 fdeja
= super(ComplexHermitianEJA
, self
)
1196 return fdeja
.__init
__(field
,
1202 def inner_product(self
, x
, y
):
1203 # Since a+bi on the diagonal is represented as
1208 # we'll double-count the "a" entries if we take the trace of
1210 return _matrix_ip(x
,y
)/2
1213 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1215 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1216 matrices, the usual symmetric Jordan product, and the
1217 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1222 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1226 The dimension of this algebra is `n^2`::
1228 sage: set_random_seed()
1229 sage: n = ZZ.random_element(1,5)
1230 sage: J = QuaternionHermitianEJA(n)
1231 sage: J.dimension() == 2*(n^2) - n
1234 The Jordan multiplication is what we think it is::
1236 sage: set_random_seed()
1237 sage: n = ZZ.random_element(1,5)
1238 sage: J = QuaternionHermitianEJA(n)
1239 sage: x = J.random_element()
1240 sage: y = J.random_element()
1241 sage: actual = (x*y).natural_representation()
1242 sage: X = x.natural_representation()
1243 sage: Y = y.natural_representation()
1244 sage: expected = (X*Y + Y*X)/2
1245 sage: actual == expected
1247 sage: J(expected) == x*y
1251 def __init__(self
, n
, field
=QQ
):
1252 S
= _quaternion_hermitian_basis(n
)
1253 Qs
= _multiplication_table_from_matrix_basis(S
)
1255 fdeja
= super(QuaternionHermitianEJA
, self
)
1256 return fdeja
.__init
__(field
,
1261 def inner_product(self
, x
, y
):
1262 # Since a+bi+cj+dk on the diagonal is represented as
1264 # a + bi +cj + dk = [ a b c d]
1269 # we'll quadruple-count the "a" entries if we take the trace of
1271 return _matrix_ip(x
,y
)/4
1274 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1276 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1277 with the usual inner product and jordan product ``x*y =
1278 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1283 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1287 This multiplication table can be verified by hand::
1289 sage: J = JordanSpinEJA(4)
1290 sage: e0,e1,e2,e3 = J.gens()
1307 def __init__(self
, n
, field
=QQ
):
1308 V
= VectorSpace(field
, n
)
1309 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1319 z0
= x
.inner_product(y
)
1320 zbar
= y0
*xbar
+ x0
*ybar
1321 z
= V([z0
] + zbar
.list())
1322 mult_table
[i
][j
] = z
1324 # The rank of the spin algebra is two, unless we're in a
1325 # one-dimensional ambient space (because the rank is bounded by
1326 # the ambient dimension).
1327 fdeja
= super(JordanSpinEJA
, self
)
1328 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2))
1330 def inner_product(self
, x
, y
):
1331 return _usual_ip(x
,y
)