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 # Turn my vector space into a module so that "vectors" can
263 # have multivatiate polynomial entries.
264 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
265 R
= PolynomialRing(self
.base_ring(), names
)
266 V
= self
.vector_space().change_ring(R
)
268 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
272 # And figure out the "left multiplication by x" matrix in
275 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
276 for i
in range(n
) ] # don't recompute these!
278 ek
= self
.monomial(k
).to_vector()
280 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
281 for i
in range(n
) ) )
282 Lx
= matrix
.column(R
, lmbx_cols
)
284 # Now we can compute powers of x "symbolically"
285 x_powers
= [self
.one().to_vector(), x
]
286 for d
in range(2, r
+1):
287 x_powers
.append( Lx
*(x_powers
[-1]) )
289 idmat
= matrix
.identity(R
, n
)
291 W
= self
._charpoly
_basis
_space
()
292 W
= W
.change_ring(R
.fraction_field())
294 # Starting with the standard coordinates x = (X1,X2,...,Xn)
295 # and then converting the entries to W-coordinates allows us
296 # to pass in the standard coordinates to the charpoly and get
297 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
300 # W.coordinates(x^2) eval'd at (standard z-coords)
304 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
306 # We want the middle equivalent thing in our matrix, but use
307 # the first equivalent thing instead so that we can pass in
308 # standard coordinates.
309 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
310 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
311 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
312 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
316 def characteristic_polynomial(self
):
318 Return a characteristic polynomial that works for all elements
321 The resulting polynomial has `n+1` variables, where `n` is the
322 dimension of this algebra. The first `n` variables correspond to
323 the coordinates of an algebra element: when evaluated at the
324 coordinates of an algebra element with respect to a certain
325 basis, the result is a univariate polynomial (in the one
326 remaining variable ``t``), namely the characteristic polynomial
331 sage: from mjo.eja.eja_algebra import JordanSpinEJA
335 The characteristic polynomial in the spin algebra is given in
336 Alizadeh, Example 11.11::
338 sage: J = JordanSpinEJA(3)
339 sage: p = J.characteristic_polynomial(); p
340 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
341 sage: xvec = J.one().to_vector()
349 # The list of coefficient polynomials a_1, a_2, ..., a_n.
350 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
352 # We go to a bit of trouble here to reorder the
353 # indeterminates, so that it's easier to evaluate the
354 # characteristic polynomial at x's coordinates and get back
355 # something in terms of t, which is what we want.
357 S
= PolynomialRing(self
.base_ring(),'t')
359 S
= PolynomialRing(S
, R
.variable_names())
362 # Note: all entries past the rth should be zero. The
363 # coefficient of the highest power (x^r) is 1, but it doesn't
364 # appear in the solution vector which contains coefficients
365 # for the other powers (to make them sum to x^r).
367 a
[r
] = 1 # corresponds to x^r
369 # When the rank is equal to the dimension, trying to
370 # assign a[r] goes out-of-bounds.
371 a
.append(1) # corresponds to x^r
373 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
376 def inner_product(self
, x
, y
):
378 The inner product associated with this Euclidean Jordan algebra.
380 Defaults to the trace inner product, but can be overridden by
381 subclasses if they are sure that the necessary properties are
386 sage: from mjo.eja.eja_algebra import random_eja
390 The inner product must satisfy its axiom for this algebra to truly
391 be a Euclidean Jordan Algebra::
393 sage: set_random_seed()
394 sage: J = random_eja()
395 sage: x = J.random_element()
396 sage: y = J.random_element()
397 sage: z = J.random_element()
398 sage: (x*y).inner_product(z) == y.inner_product(x*z)
402 if (not x
in self
) or (not y
in self
):
403 raise TypeError("arguments must live in this algebra")
404 return x
.trace_inner_product(y
)
407 def multiplication_table(self
):
409 Return a visual representation of this algebra's multiplication
410 table (on basis elements).
414 sage: from mjo.eja.eja_algebra import JordanSpinEJA
418 sage: J = JordanSpinEJA(4)
419 sage: J.multiplication_table()
420 +----++----+----+----+----+
421 | * || e0 | e1 | e2 | e3 |
422 +====++====+====+====+====+
423 | e0 || e0 | e1 | e2 | e3 |
424 +----++----+----+----+----+
425 | e1 || e1 | e0 | 0 | 0 |
426 +----++----+----+----+----+
427 | e2 || e2 | 0 | e0 | 0 |
428 +----++----+----+----+----+
429 | e3 || e3 | 0 | 0 | e0 |
430 +----++----+----+----+----+
433 M
= list(self
._multiplication
_table
) # copy
434 for i
in range(len(M
)):
435 # M had better be "square"
436 M
[i
] = [self
.monomial(i
)] + M
[i
]
437 M
= [["*"] + list(self
.gens())] + M
438 return table(M
, header_row
=True, header_column
=True, frame
=True)
441 def natural_basis(self
):
443 Return a more-natural representation of this algebra's basis.
445 Every finite-dimensional Euclidean Jordan Algebra is a direct
446 sum of five simple algebras, four of which comprise Hermitian
447 matrices. This method returns the original "natural" basis
448 for our underlying vector space. (Typically, the natural basis
449 is used to construct the multiplication table in the first place.)
451 Note that this will always return a matrix. The standard basis
452 in `R^n` will be returned as `n`-by-`1` column matrices.
456 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
457 ....: RealSymmetricEJA)
461 sage: J = RealSymmetricEJA(2)
463 Finite family {0: e0, 1: e1, 2: e2}
464 sage: J.natural_basis()
472 sage: J = JordanSpinEJA(2)
474 Finite family {0: e0, 1: e1}
475 sage: J.natural_basis()
482 if self
._natural
_basis
is None:
483 return tuple( b
.to_vector().column() for b
in self
.basis() )
485 return self
._natural
_basis
491 Return the unit element of this algebra.
495 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
500 sage: J = RealCartesianProductEJA(5)
502 e0 + e1 + e2 + e3 + e4
506 The identity element acts like the identity::
508 sage: set_random_seed()
509 sage: J = random_eja()
510 sage: x = J.random_element()
511 sage: J.one()*x == x and x*J.one() == x
514 The matrix of the unit element's operator is the identity::
516 sage: set_random_seed()
517 sage: J = random_eja()
518 sage: actual = J.one().operator().matrix()
519 sage: expected = matrix.identity(J.base_ring(), J.dimension())
520 sage: actual == expected
524 # We can brute-force compute the matrices of the operators
525 # that correspond to the basis elements of this algebra.
526 # If some linear combination of those basis elements is the
527 # algebra identity, then the same linear combination of
528 # their matrices has to be the identity matrix.
530 # Of course, matrices aren't vectors in sage, so we have to
531 # appeal to the "long vectors" isometry.
532 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
534 # Now we use basis linear algebra to find the coefficients,
535 # of the matrices-as-vectors-linear-combination, which should
536 # work for the original algebra basis too.
537 A
= matrix
.column(self
.base_ring(), oper_vecs
)
539 # We used the isometry on the left-hand side already, but we
540 # still need to do it for the right-hand side. Recall that we
541 # wanted something that summed to the identity matrix.
542 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
544 # Now if there's an identity element in the algebra, this should work.
545 coeffs
= A
.solve_right(b
)
546 return self
.linear_combination(zip(self
.gens(), coeffs
))
551 Return the rank of this EJA.
555 The author knows of no algorithm to compute the rank of an EJA
556 where only the multiplication table is known. In lieu of one, we
557 require the rank to be specified when the algebra is created,
558 and simply pass along that number here.
562 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
563 ....: RealSymmetricEJA,
564 ....: ComplexHermitianEJA,
565 ....: QuaternionHermitianEJA,
570 The rank of the Jordan spin algebra is always two::
572 sage: JordanSpinEJA(2).rank()
574 sage: JordanSpinEJA(3).rank()
576 sage: JordanSpinEJA(4).rank()
579 The rank of the `n`-by-`n` Hermitian real, complex, or
580 quaternion matrices is `n`::
582 sage: RealSymmetricEJA(2).rank()
584 sage: ComplexHermitianEJA(2).rank()
586 sage: QuaternionHermitianEJA(2).rank()
588 sage: RealSymmetricEJA(5).rank()
590 sage: ComplexHermitianEJA(5).rank()
592 sage: QuaternionHermitianEJA(5).rank()
597 Ensure that every EJA that we know how to construct has a
598 positive integer rank::
600 sage: set_random_seed()
601 sage: r = random_eja().rank()
602 sage: r in ZZ and r > 0
609 def vector_space(self
):
611 Return the vector space that underlies this algebra.
615 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
619 sage: J = RealSymmetricEJA(2)
620 sage: J.vector_space()
621 Vector space of dimension 3 over Rational Field
624 return self
.zero().to_vector().parent().ambient_vector_space()
627 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
630 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
632 Return the Euclidean Jordan Algebra corresponding to the set
633 `R^n` under the Hadamard product.
635 Note: this is nothing more than the Cartesian product of ``n``
636 copies of the spin algebra. Once Cartesian product algebras
637 are implemented, this can go.
641 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
645 This multiplication table can be verified by hand::
647 sage: J = RealCartesianProductEJA(3)
648 sage: e0,e1,e2 = J.gens()
664 We can change the generator prefix::
666 sage: RealCartesianProductEJA(3, prefix='r').gens()
670 def __init__(self
, n
, field
=QQ
, **kwargs
):
671 V
= VectorSpace(field
, n
)
672 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
675 fdeja
= super(RealCartesianProductEJA
, self
)
676 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
678 def inner_product(self
, x
, y
):
679 return _usual_ip(x
,y
)
684 Return a "random" finite-dimensional Euclidean Jordan Algebra.
688 For now, we choose a random natural number ``n`` (greater than zero)
689 and then give you back one of the following:
691 * The cartesian product of the rational numbers ``n`` times; this is
692 ``QQ^n`` with the Hadamard product.
694 * The Jordan spin algebra on ``QQ^n``.
696 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
699 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
700 in the space of ``2n``-by-``2n`` real symmetric matrices.
702 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
703 in the space of ``4n``-by-``4n`` real symmetric matrices.
705 Later this might be extended to return Cartesian products of the
710 sage: from mjo.eja.eja_algebra import random_eja
715 Euclidean Jordan algebra of dimension...
719 # The max_n component lets us choose different upper bounds on the
720 # value "n" that gets passed to the constructor. This is needed
721 # because e.g. R^{10} is reasonable to test, while the Hermitian
722 # 10-by-10 quaternion matrices are not.
723 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
725 (RealSymmetricEJA
, 5),
726 (ComplexHermitianEJA
, 4),
727 (QuaternionHermitianEJA
, 3)])
728 n
= ZZ
.random_element(1, max_n
)
729 return constructor(n
, field
=QQ
)
733 def _real_symmetric_basis(n
, field
=QQ
):
735 Return a basis for the space of real symmetric n-by-n matrices.
737 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
741 for j
in xrange(i
+1):
742 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
746 # Beware, orthogonal but not normalized!
747 Sij
= Eij
+ Eij
.transpose()
752 def _complex_hermitian_basis(n
, field
=QQ
):
754 Returns a basis for the space of complex Hermitian n-by-n matrices.
758 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
762 sage: set_random_seed()
763 sage: n = ZZ.random_element(1,5)
764 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
768 F
= QuadraticField(-1, 'I')
771 # This is like the symmetric case, but we need to be careful:
773 # * We want conjugate-symmetry, not just symmetry.
774 # * The diagonal will (as a result) be real.
778 for j
in xrange(i
+1):
779 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
781 Sij
= _embed_complex_matrix(Eij
)
784 # Beware, orthogonal but not normalized! The second one
785 # has a minus because it's conjugated.
786 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
788 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
793 def _quaternion_hermitian_basis(n
, field
=QQ
):
795 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
799 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
803 sage: set_random_seed()
804 sage: n = ZZ.random_element(1,5)
805 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
809 Q
= QuaternionAlgebra(QQ
,-1,-1)
812 # This is like the symmetric case, but we need to be careful:
814 # * We want conjugate-symmetry, not just symmetry.
815 # * The diagonal will (as a result) be real.
819 for j
in xrange(i
+1):
820 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
822 Sij
= _embed_quaternion_matrix(Eij
)
825 # Beware, orthogonal but not normalized! The second,
826 # third, and fourth ones have a minus because they're
828 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
830 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
832 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
834 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
840 def _multiplication_table_from_matrix_basis(basis
):
842 At least three of the five simple Euclidean Jordan algebras have the
843 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
844 multiplication on the right is matrix multiplication. Given a basis
845 for the underlying matrix space, this function returns a
846 multiplication table (obtained by looping through the basis
847 elements) for an algebra of those matrices.
849 # In S^2, for example, we nominally have four coordinates even
850 # though the space is of dimension three only. The vector space V
851 # is supposed to hold the entire long vector, and the subspace W
852 # of V will be spanned by the vectors that arise from symmetric
853 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
854 field
= basis
[0].base_ring()
855 dimension
= basis
[0].nrows()
857 V
= VectorSpace(field
, dimension
**2)
858 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
860 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
863 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
864 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
869 def _embed_complex_matrix(M
):
871 Embed the n-by-n complex matrix ``M`` into the space of real
872 matrices of size 2n-by-2n via the map the sends each entry `z = a +
873 bi` to the block matrix ``[[a,b],[-b,a]]``.
877 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
881 sage: F = QuadraticField(-1,'i')
882 sage: x1 = F(4 - 2*i)
883 sage: x2 = F(1 + 2*i)
886 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
887 sage: _embed_complex_matrix(M)
896 Embedding is a homomorphism (isomorphism, in fact)::
898 sage: set_random_seed()
899 sage: n = ZZ.random_element(5)
900 sage: F = QuadraticField(-1, 'i')
901 sage: X = random_matrix(F, n)
902 sage: Y = random_matrix(F, n)
903 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
904 sage: expected = _embed_complex_matrix(X*Y)
905 sage: actual == expected
911 raise ValueError("the matrix 'M' must be square")
912 field
= M
.base_ring()
917 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
919 # We can drop the imaginaries here.
920 return matrix
.block(field
.base_ring(), n
, blocks
)
923 def _unembed_complex_matrix(M
):
925 The inverse of _embed_complex_matrix().
929 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
930 ....: _unembed_complex_matrix)
934 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
935 ....: [-2, 1, -4, 3],
936 ....: [ 9, 10, 11, 12],
937 ....: [-10, 9, -12, 11] ])
938 sage: _unembed_complex_matrix(A)
940 [ 10*i + 9 12*i + 11]
944 Unembedding is the inverse of embedding::
946 sage: set_random_seed()
947 sage: F = QuadraticField(-1, 'i')
948 sage: M = random_matrix(F, 3)
949 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
955 raise ValueError("the matrix 'M' must be square")
956 if not n
.mod(2).is_zero():
957 raise ValueError("the matrix 'M' must be a complex embedding")
959 F
= QuadraticField(-1, 'i')
962 # Go top-left to bottom-right (reading order), converting every
963 # 2-by-2 block we see to a single complex element.
965 for k
in xrange(n
/2):
966 for j
in xrange(n
/2):
967 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
968 if submat
[0,0] != submat
[1,1]:
969 raise ValueError('bad on-diagonal submatrix')
970 if submat
[0,1] != -submat
[1,0]:
971 raise ValueError('bad off-diagonal submatrix')
972 z
= submat
[0,0] + submat
[0,1]*i
975 return matrix(F
, n
/2, elements
)
978 def _embed_quaternion_matrix(M
):
980 Embed the n-by-n quaternion matrix ``M`` into the space of real
981 matrices of size 4n-by-4n by first sending each quaternion entry
982 `z = a + bi + cj + dk` to the block-complex matrix
983 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
988 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
992 sage: Q = QuaternionAlgebra(QQ,-1,-1)
993 sage: i,j,k = Q.gens()
994 sage: x = 1 + 2*i + 3*j + 4*k
995 sage: M = matrix(Q, 1, [[x]])
996 sage: _embed_quaternion_matrix(M)
1002 Embedding is a homomorphism (isomorphism, in fact)::
1004 sage: set_random_seed()
1005 sage: n = ZZ.random_element(5)
1006 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1007 sage: X = random_matrix(Q, n)
1008 sage: Y = random_matrix(Q, n)
1009 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1010 sage: expected = _embed_quaternion_matrix(X*Y)
1011 sage: actual == expected
1015 quaternions
= M
.base_ring()
1018 raise ValueError("the matrix 'M' must be square")
1020 F
= QuadraticField(-1, 'i')
1025 t
= z
.coefficient_tuple()
1030 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1031 [-c
+ d
*i
, a
- b
*i
]])
1032 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1034 # We should have real entries by now, so use the realest field
1035 # we've got for the return value.
1036 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1039 def _unembed_quaternion_matrix(M
):
1041 The inverse of _embed_quaternion_matrix().
1045 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1046 ....: _unembed_quaternion_matrix)
1050 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1051 ....: [-2, 1, -4, 3],
1052 ....: [-3, 4, 1, -2],
1053 ....: [-4, -3, 2, 1]])
1054 sage: _unembed_quaternion_matrix(M)
1055 [1 + 2*i + 3*j + 4*k]
1059 Unembedding is the inverse of embedding::
1061 sage: set_random_seed()
1062 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1063 sage: M = random_matrix(Q, 3)
1064 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1070 raise ValueError("the matrix 'M' must be square")
1071 if not n
.mod(4).is_zero():
1072 raise ValueError("the matrix 'M' must be a complex embedding")
1074 Q
= QuaternionAlgebra(QQ
,-1,-1)
1077 # Go top-left to bottom-right (reading order), converting every
1078 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1081 for l
in xrange(n
/4):
1082 for m
in xrange(n
/4):
1083 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1084 if submat
[0,0] != submat
[1,1].conjugate():
1085 raise ValueError('bad on-diagonal submatrix')
1086 if submat
[0,1] != -submat
[1,0].conjugate():
1087 raise ValueError('bad off-diagonal submatrix')
1088 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1089 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1092 return matrix(Q
, n
/4, elements
)
1095 # The usual inner product on R^n.
1097 return x
.to_vector().inner_product(y
.to_vector())
1099 # The inner product used for the real symmetric simple EJA.
1100 # We keep it as a separate function because e.g. the complex
1101 # algebra uses the same inner product, except divided by 2.
1102 def _matrix_ip(X
,Y
):
1103 X_mat
= X
.natural_representation()
1104 Y_mat
= Y
.natural_representation()
1105 return (X_mat
*Y_mat
).trace()
1108 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1110 The rank-n simple EJA consisting of real symmetric n-by-n
1111 matrices, the usual symmetric Jordan product, and the trace inner
1112 product. It has dimension `(n^2 + n)/2` over the reals.
1116 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1120 sage: J = RealSymmetricEJA(2)
1121 sage: e0, e1, e2 = J.gens()
1131 The dimension of this algebra is `(n^2 + n) / 2`::
1133 sage: set_random_seed()
1134 sage: n = ZZ.random_element(1,5)
1135 sage: J = RealSymmetricEJA(n)
1136 sage: J.dimension() == (n^2 + n)/2
1139 The Jordan multiplication is what we think it is::
1141 sage: set_random_seed()
1142 sage: n = ZZ.random_element(1,5)
1143 sage: J = RealSymmetricEJA(n)
1144 sage: x = J.random_element()
1145 sage: y = J.random_element()
1146 sage: actual = (x*y).natural_representation()
1147 sage: X = x.natural_representation()
1148 sage: Y = y.natural_representation()
1149 sage: expected = (X*Y + Y*X)/2
1150 sage: actual == expected
1152 sage: J(expected) == x*y
1155 We can change the generator prefix::
1157 sage: RealSymmetricEJA(3, prefix='q').gens()
1158 (q0, q1, q2, q3, q4, q5)
1161 def __init__(self
, n
, field
=QQ
, **kwargs
):
1162 S
= _real_symmetric_basis(n
, field
=field
)
1163 Qs
= _multiplication_table_from_matrix_basis(S
)
1165 fdeja
= super(RealSymmetricEJA
, self
)
1166 return fdeja
.__init
__(field
,
1172 def inner_product(self
, x
, y
):
1173 return _matrix_ip(x
,y
)
1176 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1178 The rank-n simple EJA consisting of complex Hermitian n-by-n
1179 matrices over the real numbers, the usual symmetric Jordan product,
1180 and the real-part-of-trace inner product. It has dimension `n^2` over
1185 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1189 The dimension of this algebra is `n^2`::
1191 sage: set_random_seed()
1192 sage: n = ZZ.random_element(1,5)
1193 sage: J = ComplexHermitianEJA(n)
1194 sage: J.dimension() == n^2
1197 The Jordan multiplication is what we think it is::
1199 sage: set_random_seed()
1200 sage: n = ZZ.random_element(1,5)
1201 sage: J = ComplexHermitianEJA(n)
1202 sage: x = J.random_element()
1203 sage: y = J.random_element()
1204 sage: actual = (x*y).natural_representation()
1205 sage: X = x.natural_representation()
1206 sage: Y = y.natural_representation()
1207 sage: expected = (X*Y + Y*X)/2
1208 sage: actual == expected
1210 sage: J(expected) == x*y
1213 We can change the generator prefix::
1215 sage: ComplexHermitianEJA(2, prefix='z').gens()
1219 def __init__(self
, n
, field
=QQ
, **kwargs
):
1220 S
= _complex_hermitian_basis(n
)
1221 Qs
= _multiplication_table_from_matrix_basis(S
)
1223 fdeja
= super(ComplexHermitianEJA
, self
)
1224 return fdeja
.__init
__(field
,
1231 def inner_product(self
, x
, y
):
1232 # Since a+bi on the diagonal is represented as
1237 # we'll double-count the "a" entries if we take the trace of
1239 return _matrix_ip(x
,y
)/2
1242 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1244 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1245 matrices, the usual symmetric Jordan product, and the
1246 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1251 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1255 The dimension of this algebra is `n^2`::
1257 sage: set_random_seed()
1258 sage: n = ZZ.random_element(1,5)
1259 sage: J = QuaternionHermitianEJA(n)
1260 sage: J.dimension() == 2*(n^2) - n
1263 The Jordan multiplication is what we think it is::
1265 sage: set_random_seed()
1266 sage: n = ZZ.random_element(1,5)
1267 sage: J = QuaternionHermitianEJA(n)
1268 sage: x = J.random_element()
1269 sage: y = J.random_element()
1270 sage: actual = (x*y).natural_representation()
1271 sage: X = x.natural_representation()
1272 sage: Y = y.natural_representation()
1273 sage: expected = (X*Y + Y*X)/2
1274 sage: actual == expected
1276 sage: J(expected) == x*y
1279 We can change the generator prefix::
1281 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1282 (a0, a1, a2, a3, a4, a5)
1285 def __init__(self
, n
, field
=QQ
, **kwargs
):
1286 S
= _quaternion_hermitian_basis(n
)
1287 Qs
= _multiplication_table_from_matrix_basis(S
)
1289 fdeja
= super(QuaternionHermitianEJA
, self
)
1290 return fdeja
.__init
__(field
,
1296 def inner_product(self
, x
, y
):
1297 # Since a+bi+cj+dk on the diagonal is represented as
1299 # a + bi +cj + dk = [ a b c d]
1304 # we'll quadruple-count the "a" entries if we take the trace of
1306 return _matrix_ip(x
,y
)/4
1309 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1311 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1312 with the usual inner product and jordan product ``x*y =
1313 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1318 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1322 This multiplication table can be verified by hand::
1324 sage: J = JordanSpinEJA(4)
1325 sage: e0,e1,e2,e3 = J.gens()
1341 We can change the generator prefix::
1343 sage: JordanSpinEJA(2, prefix='B').gens()
1347 def __init__(self
, n
, field
=QQ
, **kwargs
):
1348 V
= VectorSpace(field
, n
)
1349 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1359 z0
= x
.inner_product(y
)
1360 zbar
= y0
*xbar
+ x0
*ybar
1361 z
= V([z0
] + zbar
.list())
1362 mult_table
[i
][j
] = z
1364 # The rank of the spin algebra is two, unless we're in a
1365 # one-dimensional ambient space (because the rank is bounded by
1366 # the ambient dimension).
1367 fdeja
= super(JordanSpinEJA
, self
)
1368 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1370 def inner_product(self
, x
, y
):
1371 return _usual_ip(x
,y
)