2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
9 from sage
.categories
.finite_dimensional_algebras_with_basis
import FiniteDimensionalAlgebrasWithBasis
10 from sage
.combinat
.free_module
import CombinatorialFreeModule
11 from sage
.matrix
.constructor
import matrix
12 from sage
.misc
.cachefunc
import cached_method
13 from sage
.misc
.prandom
import choice
14 from sage
.modules
.free_module
import VectorSpace
15 from sage
.rings
.integer_ring
import ZZ
16 from sage
.rings
.number_field
.number_field
import QuadraticField
17 from sage
.rings
.polynomial
.polynomial_ring_constructor
import PolynomialRing
18 from sage
.rings
.rational_field
import QQ
19 from sage
.structure
.element
import is_Matrix
21 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
22 from mjo
.eja
.eja_utils
import _mat2vec
24 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
35 sage: from mjo.eja.eja_algebra import random_eja
39 By definition, Jordan multiplication commutes::
41 sage: set_random_seed()
42 sage: J = random_eja()
43 sage: x = J.random_element()
44 sage: y = J.random_element()
50 self
._natural
_basis
= natural_basis
53 category
= FiniteDimensionalAlgebrasWithBasis(field
).Unital()
54 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
56 range(len(mult_table
)),
59 self
.print_options(bracket
='')
61 # The multiplication table we're given is necessarily in terms
62 # of vectors, because we don't have an algebra yet for
63 # anything to be an element of. However, it's faster in the
64 # long run to have the multiplication table be in terms of
65 # algebra elements. We do this after calling the superclass
66 # constructor so that from_vector() knows what to do.
67 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
68 for ls
in mult_table
]
71 def _element_constructor_(self
, elt
):
73 Construct an element of this algebra from its natural
76 This gets called only after the parent element _call_ method
77 fails to find a coercion for the argument.
81 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
82 ....: RealCartesianProductEJA,
83 ....: RealSymmetricEJA)
87 The identity in `S^n` is converted to the identity in the EJA::
89 sage: J = RealSymmetricEJA(3)
90 sage: I = matrix.identity(QQ,3)
94 This skew-symmetric matrix can't be represented in the EJA::
96 sage: J = RealSymmetricEJA(3)
97 sage: A = matrix(QQ,3, lambda i,j: i-j)
99 Traceback (most recent call last):
101 ArithmeticError: vector is not in free module
105 Ensure that we can convert any element of the two non-matrix
106 simple algebras (whose natural representations are their usual
107 vector representations) back and forth faithfully::
109 sage: set_random_seed()
110 sage: J = RealCartesianProductEJA(5)
111 sage: x = J.random_element()
112 sage: J(x.to_vector().column()) == x
114 sage: J = JordanSpinEJA(5)
115 sage: x = J.random_element()
116 sage: J(x.to_vector().column()) == x
120 natural_basis
= self
.natural_basis()
121 if elt
not in natural_basis
[0].matrix_space():
122 raise ValueError("not a naturally-represented algebra element")
124 # Thanks for nothing! Matrix spaces aren't vector
125 # spaces in Sage, so we have to figure out its
126 # natural-basis coordinates ourselves.
127 V
= VectorSpace(elt
.base_ring(), elt
.nrows()*elt
.ncols())
128 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
129 coords
= W
.coordinate_vector(_mat2vec(elt
))
130 return self
.from_vector(coords
)
135 Return a string representation of ``self``.
139 sage: from mjo.eja.eja_algebra import JordanSpinEJA
143 Ensure that it says what we think it says::
145 sage: JordanSpinEJA(2, field=QQ)
146 Euclidean Jordan algebra of degree 2 over Rational Field
147 sage: JordanSpinEJA(3, field=RDF)
148 Euclidean Jordan algebra of degree 3 over Real Double Field
151 # TODO: change this to say "dimension" and fix all the tests.
152 fmt
= "Euclidean Jordan algebra of degree {} over {}"
153 return fmt
.format(self
.dimension(), self
.base_ring())
155 def product_on_basis(self
, i
, j
):
156 return self
._multiplication
_table
[i
][j
]
158 def _a_regular_element(self
):
160 Guess a regular element. Needed to compute the basis for our
161 characteristic polynomial coefficients.
165 sage: from mjo.eja.eja_algebra import random_eja
169 Ensure that this hacky method succeeds for every algebra that we
170 know how to construct::
172 sage: set_random_seed()
173 sage: J = random_eja()
174 sage: J._a_regular_element().is_regular()
179 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
180 if not z
.is_regular():
181 raise ValueError("don't know a regular element")
186 def _charpoly_basis_space(self
):
188 Return the vector space spanned by the basis used in our
189 characteristic polynomial coefficients. This is used not only to
190 compute those coefficients, but also any time we need to
191 evaluate the coefficients (like when we compute the trace or
194 z
= self
._a
_regular
_element
()
195 V
= self
.vector_space()
196 V1
= V
.span_of_basis( (z
**k
).to_vector() for k
in range(self
.rank()) )
197 b
= (V1
.basis() + V1
.complement().basis())
198 return V
.span_of_basis(b
)
202 def _charpoly_coeff(self
, i
):
204 Return the coefficient polynomial "a_{i}" of this algebra's
205 general characteristic polynomial.
207 Having this be a separate cached method lets us compute and
208 store the trace/determinant (a_{r-1} and a_{0} respectively)
209 separate from the entire characteristic polynomial.
211 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
212 R
= A_of_x
.base_ring()
214 # Guaranteed by theory
217 # Danger: the in-place modification is done for performance
218 # reasons (reconstructing a matrix with huge polynomial
219 # entries is slow), but I don't know how cached_method works,
220 # so it's highly possible that we're modifying some global
221 # list variable by reference, here. In other words, you
222 # probably shouldn't call this method twice on the same
223 # algebra, at the same time, in two threads
224 Ai_orig
= A_of_x
.column(i
)
225 A_of_x
.set_column(i
,xr
)
226 numerator
= A_of_x
.det()
227 A_of_x
.set_column(i
,Ai_orig
)
229 # We're relying on the theory here to ensure that each a_i is
230 # indeed back in R, and the added negative signs are to make
231 # the whole charpoly expression sum to zero.
232 return R(-numerator
/detA
)
236 def _charpoly_matrix_system(self
):
238 Compute the matrix whose entries A_ij are polynomials in
239 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
240 corresponding to `x^r` and the determinent of the matrix A =
241 [A_ij]. In other words, all of the fixed (cachable) data needed
242 to compute the coefficients of the characteristic polynomial.
247 # Construct a new algebra over a multivariate polynomial ring...
248 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
249 R
= PolynomialRing(self
.base_ring(), names
)
250 # Hack around the fact that our multiplication table is in terms of
251 # algebra elements but the constructor wants it in terms of vectors.
252 vmt
= [ tuple(map(lambda x
: x
.to_vector(), ls
))
253 for ls
in self
._multiplication
_table
]
254 J
= FiniteDimensionalEuclideanJordanAlgebra(R
, tuple(vmt
), r
)
256 idmat
= matrix
.identity(J
.base_ring(), n
)
258 W
= self
._charpoly
_basis
_space
()
259 W
= W
.change_ring(R
.fraction_field())
261 # Starting with the standard coordinates x = (X1,X2,...,Xn)
262 # and then converting the entries to W-coordinates allows us
263 # to pass in the standard coordinates to the charpoly and get
264 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
267 # W.coordinates(x^2) eval'd at (standard z-coords)
271 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
273 # We want the middle equivalent thing in our matrix, but use
274 # the first equivalent thing instead so that we can pass in
275 # standard coordinates.
276 x
= J
.from_vector(W(R
.gens()))
278 # Handle the zeroth power separately, because computing
279 # the unit element in J is mathematically suspect.
280 x0
= W
.coordinate_vector(self
.one().to_vector())
282 l1
+= [ W
.coordinate_vector((x
**k
).to_vector()).column()
283 for k
in range(1,r
) ]
284 l2
= [idmat
.column(k
-1).column() for k
in range(r
+1, n
+1)]
285 A_of_x
= matrix
.block(R
, 1, n
, (l1
+ l2
))
286 xr
= W
.coordinate_vector((x
**r
).to_vector())
287 return (A_of_x
, x
, xr
, A_of_x
.det())
291 def characteristic_polynomial(self
):
293 Return a characteristic polynomial that works for all elements
296 The resulting polynomial has `n+1` variables, where `n` is the
297 dimension of this algebra. The first `n` variables correspond to
298 the coordinates of an algebra element: when evaluated at the
299 coordinates of an algebra element with respect to a certain
300 basis, the result is a univariate polynomial (in the one
301 remaining variable ``t``), namely the characteristic polynomial
306 sage: from mjo.eja.eja_algebra import JordanSpinEJA
310 The characteristic polynomial in the spin algebra is given in
311 Alizadeh, Example 11.11::
313 sage: J = JordanSpinEJA(3)
314 sage: p = J.characteristic_polynomial(); p
315 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
316 sage: xvec = J.one().to_vector()
324 # The list of coefficient polynomials a_1, a_2, ..., a_n.
325 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
327 # We go to a bit of trouble here to reorder the
328 # indeterminates, so that it's easier to evaluate the
329 # characteristic polynomial at x's coordinates and get back
330 # something in terms of t, which is what we want.
332 S
= PolynomialRing(self
.base_ring(),'t')
334 S
= PolynomialRing(S
, R
.variable_names())
337 # Note: all entries past the rth should be zero. The
338 # coefficient of the highest power (x^r) is 1, but it doesn't
339 # appear in the solution vector which contains coefficients
340 # for the other powers (to make them sum to x^r).
342 a
[r
] = 1 # corresponds to x^r
344 # When the rank is equal to the dimension, trying to
345 # assign a[r] goes out-of-bounds.
346 a
.append(1) # corresponds to x^r
348 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
351 def inner_product(self
, x
, y
):
353 The inner product associated with this Euclidean Jordan algebra.
355 Defaults to the trace inner product, but can be overridden by
356 subclasses if they are sure that the necessary properties are
361 sage: from mjo.eja.eja_algebra import random_eja
365 The inner product must satisfy its axiom for this algebra to truly
366 be a Euclidean Jordan Algebra::
368 sage: set_random_seed()
369 sage: J = random_eja()
370 sage: x = J.random_element()
371 sage: y = J.random_element()
372 sage: z = J.random_element()
373 sage: (x*y).inner_product(z) == y.inner_product(x*z)
377 if (not x
in self
) or (not y
in self
):
378 raise TypeError("arguments must live in this algebra")
379 return x
.trace_inner_product(y
)
382 def multiplication_table(self
):
384 Return a readable matrix representation of this algebra's
385 multiplication table. The (i,j)th entry in the matrix contains
386 the product of the ith basis element with the jth.
388 This is not extraordinarily useful, but it overrides a superclass
389 method that would otherwise just crash and complain about the
390 algebra being infinite.
394 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
395 ....: RealCartesianProductEJA)
399 sage: J = RealCartesianProductEJA(3)
400 sage: J.multiplication_table()
407 sage: J = JordanSpinEJA(3)
408 sage: J.multiplication_table()
414 return matrix(self
._multiplication
_table
)
417 def natural_basis(self
):
419 Return a more-natural representation of this algebra's basis.
421 Every finite-dimensional Euclidean Jordan Algebra is a direct
422 sum of five simple algebras, four of which comprise Hermitian
423 matrices. This method returns the original "natural" basis
424 for our underlying vector space. (Typically, the natural basis
425 is used to construct the multiplication table in the first place.)
427 Note that this will always return a matrix. The standard basis
428 in `R^n` will be returned as `n`-by-`1` column matrices.
432 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
433 ....: RealSymmetricEJA)
437 sage: J = RealSymmetricEJA(2)
439 Finite family {0: e0, 1: e1, 2: e2}
440 sage: J.natural_basis()
448 sage: J = JordanSpinEJA(2)
450 Finite family {0: e0, 1: e1}
451 sage: J.natural_basis()
458 if self
._natural
_basis
is None:
459 return tuple( b
.to_vector().column() for b
in self
.basis() )
461 return self
._natural
_basis
467 Return the unit element of this algebra.
471 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
476 sage: J = RealCartesianProductEJA(5)
478 e0 + e1 + e2 + e3 + e4
482 The identity element acts like the identity::
484 sage: set_random_seed()
485 sage: J = random_eja()
486 sage: x = J.random_element()
487 sage: J.one()*x == x and x*J.one() == x
490 The matrix of the unit element's operator is the identity::
492 sage: set_random_seed()
493 sage: J = random_eja()
494 sage: actual = J.one().operator().matrix()
495 sage: expected = matrix.identity(J.base_ring(), J.dimension())
496 sage: actual == expected
500 # We can brute-force compute the matrices of the operators
501 # that correspond to the basis elements of this algebra.
502 # If some linear combination of those basis elements is the
503 # algebra identity, then the same linear combination of
504 # their matrices has to be the identity matrix.
506 # Of course, matrices aren't vectors in sage, so we have to
507 # appeal to the "long vectors" isometry.
508 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
510 # Now we use basis linear algebra to find the coefficients,
511 # of the matrices-as-vectors-linear-combination, which should
512 # work for the original algebra basis too.
513 A
= matrix
.column(self
.base_ring(), oper_vecs
)
515 # We used the isometry on the left-hand side already, but we
516 # still need to do it for the right-hand side. Recall that we
517 # wanted something that summed to the identity matrix.
518 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
520 # Now if there's an identity element in the algebra, this should work.
521 coeffs
= A
.solve_right(b
)
522 return self
.linear_combination(zip(self
.gens(), coeffs
))
527 Return the rank of this EJA.
531 The author knows of no algorithm to compute the rank of an EJA
532 where only the multiplication table is known. In lieu of one, we
533 require the rank to be specified when the algebra is created,
534 and simply pass along that number here.
538 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
539 ....: RealSymmetricEJA,
540 ....: ComplexHermitianEJA,
541 ....: QuaternionHermitianEJA,
546 The rank of the Jordan spin algebra is always two::
548 sage: JordanSpinEJA(2).rank()
550 sage: JordanSpinEJA(3).rank()
552 sage: JordanSpinEJA(4).rank()
555 The rank of the `n`-by-`n` Hermitian real, complex, or
556 quaternion matrices is `n`::
558 sage: RealSymmetricEJA(2).rank()
560 sage: ComplexHermitianEJA(2).rank()
562 sage: QuaternionHermitianEJA(2).rank()
564 sage: RealSymmetricEJA(5).rank()
566 sage: ComplexHermitianEJA(5).rank()
568 sage: QuaternionHermitianEJA(5).rank()
573 Ensure that every EJA that we know how to construct has a
574 positive integer rank::
576 sage: set_random_seed()
577 sage: r = random_eja().rank()
578 sage: r in ZZ and r > 0
585 def vector_space(self
):
587 Return the vector space that underlies this algebra.
591 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
595 sage: J = RealSymmetricEJA(2)
596 sage: J.vector_space()
597 Vector space of dimension 3 over Rational Field
600 return self
.zero().to_vector().parent().ambient_vector_space()
603 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
606 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
608 Return the Euclidean Jordan Algebra corresponding to the set
609 `R^n` under the Hadamard product.
611 Note: this is nothing more than the Cartesian product of ``n``
612 copies of the spin algebra. Once Cartesian product algebras
613 are implemented, this can go.
617 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
621 This multiplication table can be verified by hand::
623 sage: J = RealCartesianProductEJA(3)
624 sage: e0,e1,e2 = J.gens()
639 def __init__(self
, n
, field
=QQ
):
640 V
= VectorSpace(field
, n
)
641 mult_table
= [ [ V
.basis()[i
]*(i
== j
) for i
in range(n
) ]
644 fdeja
= super(RealCartesianProductEJA
, self
)
645 return fdeja
.__init
__(field
, mult_table
, rank
=n
)
647 def inner_product(self
, x
, y
):
648 return _usual_ip(x
,y
)
653 Return a "random" finite-dimensional Euclidean Jordan Algebra.
657 For now, we choose a random natural number ``n`` (greater than zero)
658 and then give you back one of the following:
660 * The cartesian product of the rational numbers ``n`` times; this is
661 ``QQ^n`` with the Hadamard product.
663 * The Jordan spin algebra on ``QQ^n``.
665 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
668 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
669 in the space of ``2n``-by-``2n`` real symmetric matrices.
671 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
672 in the space of ``4n``-by-``4n`` real symmetric matrices.
674 Later this might be extended to return Cartesian products of the
679 sage: from mjo.eja.eja_algebra import random_eja
684 Euclidean Jordan algebra of degree...
688 # The max_n component lets us choose different upper bounds on the
689 # value "n" that gets passed to the constructor. This is needed
690 # because e.g. R^{10} is reasonable to test, while the Hermitian
691 # 10-by-10 quaternion matrices are not.
692 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
694 (RealSymmetricEJA
, 5),
695 (ComplexHermitianEJA
, 4),
696 (QuaternionHermitianEJA
, 3)])
697 n
= ZZ
.random_element(1, max_n
)
698 return constructor(n
, field
=QQ
)
702 def _real_symmetric_basis(n
, field
=QQ
):
704 Return a basis for the space of real symmetric n-by-n matrices.
706 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
710 for j
in xrange(i
+1):
711 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
715 # Beware, orthogonal but not normalized!
716 Sij
= Eij
+ Eij
.transpose()
721 def _complex_hermitian_basis(n
, field
=QQ
):
723 Returns a basis for the space of complex Hermitian n-by-n matrices.
727 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
731 sage: set_random_seed()
732 sage: n = ZZ.random_element(1,5)
733 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
737 F
= QuadraticField(-1, 'I')
740 # This is like the symmetric case, but we need to be careful:
742 # * We want conjugate-symmetry, not just symmetry.
743 # * The diagonal will (as a result) be real.
747 for j
in xrange(i
+1):
748 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
750 Sij
= _embed_complex_matrix(Eij
)
753 # Beware, orthogonal but not normalized! The second one
754 # has a minus because it's conjugated.
755 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
757 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
762 def _quaternion_hermitian_basis(n
, field
=QQ
):
764 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
768 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
772 sage: set_random_seed()
773 sage: n = ZZ.random_element(1,5)
774 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
778 Q
= QuaternionAlgebra(QQ
,-1,-1)
781 # This is like the symmetric case, but we need to be careful:
783 # * We want conjugate-symmetry, not just symmetry.
784 # * The diagonal will (as a result) be real.
788 for j
in xrange(i
+1):
789 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
791 Sij
= _embed_quaternion_matrix(Eij
)
794 # Beware, orthogonal but not normalized! The second,
795 # third, and fourth ones have a minus because they're
797 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
799 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
801 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
803 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
809 def _multiplication_table_from_matrix_basis(basis
):
811 At least three of the five simple Euclidean Jordan algebras have the
812 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
813 multiplication on the right is matrix multiplication. Given a basis
814 for the underlying matrix space, this function returns a
815 multiplication table (obtained by looping through the basis
816 elements) for an algebra of those matrices.
818 # In S^2, for example, we nominally have four coordinates even
819 # though the space is of dimension three only. The vector space V
820 # is supposed to hold the entire long vector, and the subspace W
821 # of V will be spanned by the vectors that arise from symmetric
822 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
823 field
= basis
[0].base_ring()
824 dimension
= basis
[0].nrows()
826 V
= VectorSpace(field
, dimension
**2)
827 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
829 mult_table
= [[W
.zero() for i
in range(n
)] for j
in range(n
)]
832 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
833 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
838 def _embed_complex_matrix(M
):
840 Embed the n-by-n complex matrix ``M`` into the space of real
841 matrices of size 2n-by-2n via the map the sends each entry `z = a +
842 bi` to the block matrix ``[[a,b],[-b,a]]``.
846 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
850 sage: F = QuadraticField(-1,'i')
851 sage: x1 = F(4 - 2*i)
852 sage: x2 = F(1 + 2*i)
855 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
856 sage: _embed_complex_matrix(M)
865 Embedding is a homomorphism (isomorphism, in fact)::
867 sage: set_random_seed()
868 sage: n = ZZ.random_element(5)
869 sage: F = QuadraticField(-1, 'i')
870 sage: X = random_matrix(F, n)
871 sage: Y = random_matrix(F, n)
872 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
873 sage: expected = _embed_complex_matrix(X*Y)
874 sage: actual == expected
880 raise ValueError("the matrix 'M' must be square")
881 field
= M
.base_ring()
886 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
888 # We can drop the imaginaries here.
889 return matrix
.block(field
.base_ring(), n
, blocks
)
892 def _unembed_complex_matrix(M
):
894 The inverse of _embed_complex_matrix().
898 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
899 ....: _unembed_complex_matrix)
903 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
904 ....: [-2, 1, -4, 3],
905 ....: [ 9, 10, 11, 12],
906 ....: [-10, 9, -12, 11] ])
907 sage: _unembed_complex_matrix(A)
909 [ 10*i + 9 12*i + 11]
913 Unembedding is the inverse of embedding::
915 sage: set_random_seed()
916 sage: F = QuadraticField(-1, 'i')
917 sage: M = random_matrix(F, 3)
918 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
924 raise ValueError("the matrix 'M' must be square")
925 if not n
.mod(2).is_zero():
926 raise ValueError("the matrix 'M' must be a complex embedding")
928 F
= QuadraticField(-1, 'i')
931 # Go top-left to bottom-right (reading order), converting every
932 # 2-by-2 block we see to a single complex element.
934 for k
in xrange(n
/2):
935 for j
in xrange(n
/2):
936 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
937 if submat
[0,0] != submat
[1,1]:
938 raise ValueError('bad on-diagonal submatrix')
939 if submat
[0,1] != -submat
[1,0]:
940 raise ValueError('bad off-diagonal submatrix')
941 z
= submat
[0,0] + submat
[0,1]*i
944 return matrix(F
, n
/2, elements
)
947 def _embed_quaternion_matrix(M
):
949 Embed the n-by-n quaternion matrix ``M`` into the space of real
950 matrices of size 4n-by-4n by first sending each quaternion entry
951 `z = a + bi + cj + dk` to the block-complex matrix
952 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
957 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
961 sage: Q = QuaternionAlgebra(QQ,-1,-1)
962 sage: i,j,k = Q.gens()
963 sage: x = 1 + 2*i + 3*j + 4*k
964 sage: M = matrix(Q, 1, [[x]])
965 sage: _embed_quaternion_matrix(M)
971 Embedding is a homomorphism (isomorphism, in fact)::
973 sage: set_random_seed()
974 sage: n = ZZ.random_element(5)
975 sage: Q = QuaternionAlgebra(QQ,-1,-1)
976 sage: X = random_matrix(Q, n)
977 sage: Y = random_matrix(Q, n)
978 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
979 sage: expected = _embed_quaternion_matrix(X*Y)
980 sage: actual == expected
984 quaternions
= M
.base_ring()
987 raise ValueError("the matrix 'M' must be square")
989 F
= QuadraticField(-1, 'i')
994 t
= z
.coefficient_tuple()
999 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1000 [-c
+ d
*i
, a
- b
*i
]])
1001 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1003 # We should have real entries by now, so use the realest field
1004 # we've got for the return value.
1005 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1008 def _unembed_quaternion_matrix(M
):
1010 The inverse of _embed_quaternion_matrix().
1014 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1015 ....: _unembed_quaternion_matrix)
1019 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1020 ....: [-2, 1, -4, 3],
1021 ....: [-3, 4, 1, -2],
1022 ....: [-4, -3, 2, 1]])
1023 sage: _unembed_quaternion_matrix(M)
1024 [1 + 2*i + 3*j + 4*k]
1028 Unembedding is the inverse of embedding::
1030 sage: set_random_seed()
1031 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1032 sage: M = random_matrix(Q, 3)
1033 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1039 raise ValueError("the matrix 'M' must be square")
1040 if not n
.mod(4).is_zero():
1041 raise ValueError("the matrix 'M' must be a complex embedding")
1043 Q
= QuaternionAlgebra(QQ
,-1,-1)
1046 # Go top-left to bottom-right (reading order), converting every
1047 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1050 for l
in xrange(n
/4):
1051 for m
in xrange(n
/4):
1052 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1053 if submat
[0,0] != submat
[1,1].conjugate():
1054 raise ValueError('bad on-diagonal submatrix')
1055 if submat
[0,1] != -submat
[1,0].conjugate():
1056 raise ValueError('bad off-diagonal submatrix')
1057 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1058 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1061 return matrix(Q
, n
/4, elements
)
1064 # The usual inner product on R^n.
1066 return x
.to_vector().inner_product(y
.to_vector())
1068 # The inner product used for the real symmetric simple EJA.
1069 # We keep it as a separate function because e.g. the complex
1070 # algebra uses the same inner product, except divided by 2.
1071 def _matrix_ip(X
,Y
):
1072 X_mat
= X
.natural_representation()
1073 Y_mat
= Y
.natural_representation()
1074 return (X_mat
*Y_mat
).trace()
1077 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1079 The rank-n simple EJA consisting of real symmetric n-by-n
1080 matrices, the usual symmetric Jordan product, and the trace inner
1081 product. It has dimension `(n^2 + n)/2` over the reals.
1085 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1089 sage: J = RealSymmetricEJA(2)
1090 sage: e0, e1, e2 = J.gens()
1100 The dimension of this algebra is `(n^2 + n) / 2`::
1102 sage: set_random_seed()
1103 sage: n = ZZ.random_element(1,5)
1104 sage: J = RealSymmetricEJA(n)
1105 sage: J.dimension() == (n^2 + n)/2
1108 The Jordan multiplication is what we think it is::
1110 sage: set_random_seed()
1111 sage: n = ZZ.random_element(1,5)
1112 sage: J = RealSymmetricEJA(n)
1113 sage: x = J.random_element()
1114 sage: y = J.random_element()
1115 sage: actual = (x*y).natural_representation()
1116 sage: X = x.natural_representation()
1117 sage: Y = y.natural_representation()
1118 sage: expected = (X*Y + Y*X)/2
1119 sage: actual == expected
1121 sage: J(expected) == x*y
1125 def __init__(self
, n
, field
=QQ
):
1126 S
= _real_symmetric_basis(n
, field
=field
)
1127 Qs
= _multiplication_table_from_matrix_basis(S
)
1129 fdeja
= super(RealSymmetricEJA
, self
)
1130 return fdeja
.__init
__(field
,
1135 def inner_product(self
, x
, y
):
1136 return _matrix_ip(x
,y
)
1139 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1141 The rank-n simple EJA consisting of complex Hermitian n-by-n
1142 matrices over the real numbers, the usual symmetric Jordan product,
1143 and the real-part-of-trace inner product. It has dimension `n^2` over
1148 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1152 The dimension of this algebra is `n^2`::
1154 sage: set_random_seed()
1155 sage: n = ZZ.random_element(1,5)
1156 sage: J = ComplexHermitianEJA(n)
1157 sage: J.dimension() == n^2
1160 The Jordan multiplication is what we think it is::
1162 sage: set_random_seed()
1163 sage: n = ZZ.random_element(1,5)
1164 sage: J = ComplexHermitianEJA(n)
1165 sage: x = J.random_element()
1166 sage: y = J.random_element()
1167 sage: actual = (x*y).natural_representation()
1168 sage: X = x.natural_representation()
1169 sage: Y = y.natural_representation()
1170 sage: expected = (X*Y + Y*X)/2
1171 sage: actual == expected
1173 sage: J(expected) == x*y
1177 def __init__(self
, n
, field
=QQ
):
1178 S
= _complex_hermitian_basis(n
)
1179 Qs
= _multiplication_table_from_matrix_basis(S
)
1181 fdeja
= super(ComplexHermitianEJA
, self
)
1182 return fdeja
.__init
__(field
,
1188 def inner_product(self
, x
, y
):
1189 # Since a+bi on the diagonal is represented as
1194 # we'll double-count the "a" entries if we take the trace of
1196 return _matrix_ip(x
,y
)/2
1199 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1201 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1202 matrices, the usual symmetric Jordan product, and the
1203 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1208 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1212 The dimension of this algebra is `n^2`::
1214 sage: set_random_seed()
1215 sage: n = ZZ.random_element(1,5)
1216 sage: J = QuaternionHermitianEJA(n)
1217 sage: J.dimension() == 2*(n^2) - n
1220 The Jordan multiplication is what we think it is::
1222 sage: set_random_seed()
1223 sage: n = ZZ.random_element(1,5)
1224 sage: J = QuaternionHermitianEJA(n)
1225 sage: x = J.random_element()
1226 sage: y = J.random_element()
1227 sage: actual = (x*y).natural_representation()
1228 sage: X = x.natural_representation()
1229 sage: Y = y.natural_representation()
1230 sage: expected = (X*Y + Y*X)/2
1231 sage: actual == expected
1233 sage: J(expected) == x*y
1237 def __init__(self
, n
, field
=QQ
):
1238 S
= _quaternion_hermitian_basis(n
)
1239 Qs
= _multiplication_table_from_matrix_basis(S
)
1241 fdeja
= super(QuaternionHermitianEJA
, self
)
1242 return fdeja
.__init
__(field
,
1247 def inner_product(self
, x
, y
):
1248 # Since a+bi+cj+dk on the diagonal is represented as
1250 # a + bi +cj + dk = [ a b c d]
1255 # we'll quadruple-count the "a" entries if we take the trace of
1257 return _matrix_ip(x
,y
)/4
1260 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1262 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1263 with the usual inner product and jordan product ``x*y =
1264 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1269 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1273 This multiplication table can be verified by hand::
1275 sage: J = JordanSpinEJA(4)
1276 sage: e0,e1,e2,e3 = J.gens()
1293 def __init__(self
, n
, field
=QQ
):
1294 V
= VectorSpace(field
, n
)
1295 mult_table
= [[V
.zero() for i
in range(n
)] for j
in range(n
)]
1305 z0
= x
.inner_product(y
)
1306 zbar
= y0
*xbar
+ x0
*ybar
1307 z
= V([z0
] + zbar
.list())
1308 mult_table
[i
][j
] = z
1310 # The rank of the spin algebra is two, unless we're in a
1311 # one-dimensional ambient space (because the rank is bounded by
1312 # the ambient dimension).
1313 fdeja
= super(JordanSpinEJA
, self
)
1314 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2))
1316 def inner_product(self
, x
, y
):
1317 return _usual_ip(x
,y
)