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
.matrix
.matrix_space
import MatrixSpace
13 from sage
.misc
.cachefunc
import cached_method
14 from sage
.misc
.prandom
import choice
15 from sage
.misc
.table
import table
16 from sage
.modules
.free_module
import FreeModule
, VectorSpace
17 from sage
.rings
.integer_ring
import ZZ
18 from sage
.rings
.number_field
.number_field
import QuadraticField
19 from sage
.rings
.polynomial
.polynomial_ring_constructor
import PolynomialRing
20 from sage
.rings
.rational_field
import QQ
21 from sage
.structure
.element
import is_Matrix
23 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
24 from mjo
.eja
.eja_utils
import _mat2vec
26 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
27 # This is an ugly hack needed to prevent the category framework
28 # from implementing a coercion from our base ring (e.g. the
29 # rationals) into the algebra. First of all -- such a coercion is
30 # nonsense to begin with. But more importantly, it tries to do so
31 # in the category of rings, and since our algebras aren't
32 # associative they generally won't be rings.
33 _no_generic_basering_coercion
= True
45 sage: from mjo.eja.eja_algebra import random_eja
49 By definition, Jordan multiplication commutes::
51 sage: set_random_seed()
52 sage: J = random_eja()
53 sage: x = J.random_element()
54 sage: y = J.random_element()
60 self
._natural
_basis
= natural_basis
63 category
= MagmaticAlgebras(field
).FiniteDimensional()
64 category
= category
.WithBasis().Unital()
66 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
68 range(len(mult_table
)),
71 self
.print_options(bracket
='')
73 # The multiplication table we're given is necessarily in terms
74 # of vectors, because we don't have an algebra yet for
75 # anything to be an element of. However, it's faster in the
76 # long run to have the multiplication table be in terms of
77 # algebra elements. We do this after calling the superclass
78 # constructor so that from_vector() knows what to do.
79 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
80 for ls
in mult_table
]
83 def _element_constructor_(self
, elt
):
85 Construct an element of this algebra from its natural
88 This gets called only after the parent element _call_ method
89 fails to find a coercion for the argument.
93 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
94 ....: RealCartesianProductEJA,
95 ....: RealSymmetricEJA)
99 The identity in `S^n` is converted to the identity in the EJA::
101 sage: J = RealSymmetricEJA(3)
102 sage: I = matrix.identity(QQ,3)
103 sage: J(I) == J.one()
106 This skew-symmetric matrix can't be represented in the EJA::
108 sage: J = RealSymmetricEJA(3)
109 sage: A = matrix(QQ,3, lambda i,j: i-j)
111 Traceback (most recent call last):
113 ArithmeticError: vector is not in free module
117 Ensure that we can convert any element of the two non-matrix
118 simple algebras (whose natural representations are their usual
119 vector representations) back and forth faithfully::
121 sage: set_random_seed()
122 sage: J = RealCartesianProductEJA(5)
123 sage: x = J.random_element()
124 sage: J(x.to_vector().column()) == x
126 sage: J = JordanSpinEJA(5)
127 sage: x = J.random_element()
128 sage: J(x.to_vector().column()) == x
133 # The superclass implementation of random_element()
134 # needs to be able to coerce "0" into the algebra.
137 natural_basis
= self
.natural_basis()
138 if elt
not in natural_basis
[0].matrix_space():
139 raise ValueError("not a naturally-represented algebra element")
141 # Thanks for nothing! Matrix spaces aren't vector
142 # spaces in Sage, so we have to figure out its
143 # natural-basis coordinates ourselves.
144 V
= VectorSpace(elt
.base_ring(), elt
.nrows()*elt
.ncols())
145 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
146 coords
= W
.coordinate_vector(_mat2vec(elt
))
147 return self
.from_vector(coords
)
152 Return a string representation of ``self``.
156 sage: from mjo.eja.eja_algebra import JordanSpinEJA
160 Ensure that it says what we think it says::
162 sage: JordanSpinEJA(2, field=QQ)
163 Euclidean Jordan algebra of dimension 2 over Rational Field
164 sage: JordanSpinEJA(3, field=RDF)
165 Euclidean Jordan algebra of dimension 3 over Real Double Field
168 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
169 return fmt
.format(self
.dimension(), self
.base_ring())
171 def product_on_basis(self
, i
, j
):
172 return self
._multiplication
_table
[i
][j
]
174 def _a_regular_element(self
):
176 Guess a regular element. Needed to compute the basis for our
177 characteristic polynomial coefficients.
181 sage: from mjo.eja.eja_algebra import random_eja
185 Ensure that this hacky method succeeds for every algebra that we
186 know how to construct::
188 sage: set_random_seed()
189 sage: J = random_eja()
190 sage: J._a_regular_element().is_regular()
195 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
196 if not z
.is_regular():
197 raise ValueError("don't know a regular element")
202 def _charpoly_basis_space(self
):
204 Return the vector space spanned by the basis used in our
205 characteristic polynomial coefficients. This is used not only to
206 compute those coefficients, but also any time we need to
207 evaluate the coefficients (like when we compute the trace or
210 z
= self
._a
_regular
_element
()
211 # Don't use the parent vector space directly here in case this
212 # happens to be a subalgebra. In that case, we would be e.g.
213 # two-dimensional but span_of_basis() would expect three
215 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
216 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
217 V1
= V
.span_of_basis( basis
)
218 b
= (V1
.basis() + V1
.complement().basis())
219 return V
.span_of_basis(b
)
223 def _charpoly_coeff(self
, i
):
225 Return the coefficient polynomial "a_{i}" of this algebra's
226 general characteristic polynomial.
228 Having this be a separate cached method lets us compute and
229 store the trace/determinant (a_{r-1} and a_{0} respectively)
230 separate from the entire characteristic polynomial.
232 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
233 R
= A_of_x
.base_ring()
235 # Guaranteed by theory
238 # Danger: the in-place modification is done for performance
239 # reasons (reconstructing a matrix with huge polynomial
240 # entries is slow), but I don't know how cached_method works,
241 # so it's highly possible that we're modifying some global
242 # list variable by reference, here. In other words, you
243 # probably shouldn't call this method twice on the same
244 # algebra, at the same time, in two threads
245 Ai_orig
= A_of_x
.column(i
)
246 A_of_x
.set_column(i
,xr
)
247 numerator
= A_of_x
.det()
248 A_of_x
.set_column(i
,Ai_orig
)
250 # We're relying on the theory here to ensure that each a_i is
251 # indeed back in R, and the added negative signs are to make
252 # the whole charpoly expression sum to zero.
253 return R(-numerator
/detA
)
257 def _charpoly_matrix_system(self
):
259 Compute the matrix whose entries A_ij are polynomials in
260 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
261 corresponding to `x^r` and the determinent of the matrix A =
262 [A_ij]. In other words, all of the fixed (cachable) data needed
263 to compute the coefficients of the characteristic polynomial.
268 # Turn my vector space into a module so that "vectors" can
269 # have multivatiate polynomial entries.
270 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
271 R
= PolynomialRing(self
.base_ring(), names
)
273 # Using change_ring() on the parent's vector space doesn't work
274 # here because, in a subalgebra, that vector space has a basis
275 # and change_ring() tries to bring the basis along with it. And
276 # that doesn't work unless the new ring is a PID, which it usually
280 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
284 # And figure out the "left multiplication by x" matrix in
287 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
288 for i
in range(n
) ] # don't recompute these!
290 ek
= self
.monomial(k
).to_vector()
292 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
293 for i
in range(n
) ) )
294 Lx
= matrix
.column(R
, lmbx_cols
)
296 # Now we can compute powers of x "symbolically"
297 x_powers
= [self
.one().to_vector(), x
]
298 for d
in range(2, r
+1):
299 x_powers
.append( Lx
*(x_powers
[-1]) )
301 idmat
= matrix
.identity(R
, n
)
303 W
= self
._charpoly
_basis
_space
()
304 W
= W
.change_ring(R
.fraction_field())
306 # Starting with the standard coordinates x = (X1,X2,...,Xn)
307 # and then converting the entries to W-coordinates allows us
308 # to pass in the standard coordinates to the charpoly and get
309 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
312 # W.coordinates(x^2) eval'd at (standard z-coords)
316 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
318 # We want the middle equivalent thing in our matrix, but use
319 # the first equivalent thing instead so that we can pass in
320 # standard coordinates.
321 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
322 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
323 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
324 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
328 def characteristic_polynomial(self
):
330 Return a characteristic polynomial that works for all elements
333 The resulting polynomial has `n+1` variables, where `n` is the
334 dimension of this algebra. The first `n` variables correspond to
335 the coordinates of an algebra element: when evaluated at the
336 coordinates of an algebra element with respect to a certain
337 basis, the result is a univariate polynomial (in the one
338 remaining variable ``t``), namely the characteristic polynomial
343 sage: from mjo.eja.eja_algebra import JordanSpinEJA
347 The characteristic polynomial in the spin algebra is given in
348 Alizadeh, Example 11.11::
350 sage: J = JordanSpinEJA(3)
351 sage: p = J.characteristic_polynomial(); p
352 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
353 sage: xvec = J.one().to_vector()
361 # The list of coefficient polynomials a_1, a_2, ..., a_n.
362 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
364 # We go to a bit of trouble here to reorder the
365 # indeterminates, so that it's easier to evaluate the
366 # characteristic polynomial at x's coordinates and get back
367 # something in terms of t, which is what we want.
369 S
= PolynomialRing(self
.base_ring(),'t')
371 S
= PolynomialRing(S
, R
.variable_names())
374 # Note: all entries past the rth should be zero. The
375 # coefficient of the highest power (x^r) is 1, but it doesn't
376 # appear in the solution vector which contains coefficients
377 # for the other powers (to make them sum to x^r).
379 a
[r
] = 1 # corresponds to x^r
381 # When the rank is equal to the dimension, trying to
382 # assign a[r] goes out-of-bounds.
383 a
.append(1) # corresponds to x^r
385 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
388 def inner_product(self
, x
, y
):
390 The inner product associated with this Euclidean Jordan algebra.
392 Defaults to the trace inner product, but can be overridden by
393 subclasses if they are sure that the necessary properties are
398 sage: from mjo.eja.eja_algebra import random_eja
402 The inner product must satisfy its axiom for this algebra to truly
403 be a Euclidean Jordan Algebra::
405 sage: set_random_seed()
406 sage: J = random_eja()
407 sage: x = J.random_element()
408 sage: y = J.random_element()
409 sage: z = J.random_element()
410 sage: (x*y).inner_product(z) == y.inner_product(x*z)
414 if (not x
in self
) or (not y
in self
):
415 raise TypeError("arguments must live in this algebra")
416 return x
.trace_inner_product(y
)
419 def is_trivial(self
):
421 Return whether or not this algebra is trivial.
423 A trivial algebra contains only the zero element.
427 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
431 sage: J = ComplexHermitianEJA(3)
434 sage: A = J.zero().subalgebra_generated_by()
439 return self
.dimension() == 0
442 def multiplication_table(self
):
444 Return a visual representation of this algebra's multiplication
445 table (on basis elements).
449 sage: from mjo.eja.eja_algebra import JordanSpinEJA
453 sage: J = JordanSpinEJA(4)
454 sage: J.multiplication_table()
455 +----++----+----+----+----+
456 | * || e0 | e1 | e2 | e3 |
457 +====++====+====+====+====+
458 | e0 || e0 | e1 | e2 | e3 |
459 +----++----+----+----+----+
460 | e1 || e1 | e0 | 0 | 0 |
461 +----++----+----+----+----+
462 | e2 || e2 | 0 | e0 | 0 |
463 +----++----+----+----+----+
464 | e3 || e3 | 0 | 0 | e0 |
465 +----++----+----+----+----+
468 M
= list(self
._multiplication
_table
) # copy
469 for i
in range(len(M
)):
470 # M had better be "square"
471 M
[i
] = [self
.monomial(i
)] + M
[i
]
472 M
= [["*"] + list(self
.gens())] + M
473 return table(M
, header_row
=True, header_column
=True, frame
=True)
476 def natural_basis(self
):
478 Return a more-natural representation of this algebra's basis.
480 Every finite-dimensional Euclidean Jordan Algebra is a direct
481 sum of five simple algebras, four of which comprise Hermitian
482 matrices. This method returns the original "natural" basis
483 for our underlying vector space. (Typically, the natural basis
484 is used to construct the multiplication table in the first place.)
486 Note that this will always return a matrix. The standard basis
487 in `R^n` will be returned as `n`-by-`1` column matrices.
491 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
492 ....: RealSymmetricEJA)
496 sage: J = RealSymmetricEJA(2)
498 Finite family {0: e0, 1: e1, 2: e2}
499 sage: J.natural_basis()
507 sage: J = JordanSpinEJA(2)
509 Finite family {0: e0, 1: e1}
510 sage: J.natural_basis()
517 if self
._natural
_basis
is None:
518 M
= self
.natural_basis_space()
519 return tuple( M(b
.to_vector()) for b
in self
.basis() )
521 return self
._natural
_basis
524 def natural_basis_space(self
):
526 Return the matrix space in which this algebra's natural basis
529 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
530 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
532 return self
._natural
_basis
[0].matrix_space()
538 Return the unit element of this algebra.
542 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
547 sage: J = RealCartesianProductEJA(5)
549 e0 + e1 + e2 + e3 + e4
553 The identity element acts like the identity::
555 sage: set_random_seed()
556 sage: J = random_eja()
557 sage: x = J.random_element()
558 sage: J.one()*x == x and x*J.one() == x
561 The matrix of the unit element's operator is the identity::
563 sage: set_random_seed()
564 sage: J = random_eja()
565 sage: actual = J.one().operator().matrix()
566 sage: expected = matrix.identity(J.base_ring(), J.dimension())
567 sage: actual == expected
571 # We can brute-force compute the matrices of the operators
572 # that correspond to the basis elements of this algebra.
573 # If some linear combination of those basis elements is the
574 # algebra identity, then the same linear combination of
575 # their matrices has to be the identity matrix.
577 # Of course, matrices aren't vectors in sage, so we have to
578 # appeal to the "long vectors" isometry.
579 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
581 # Now we use basis linear algebra to find the coefficients,
582 # of the matrices-as-vectors-linear-combination, which should
583 # work for the original algebra basis too.
584 A
= matrix
.column(self
.base_ring(), oper_vecs
)
586 # We used the isometry on the left-hand side already, but we
587 # still need to do it for the right-hand side. Recall that we
588 # wanted something that summed to the identity matrix.
589 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
591 # Now if there's an identity element in the algebra, this should work.
592 coeffs
= A
.solve_right(b
)
593 return self
.linear_combination(zip(self
.gens(), coeffs
))
596 def random_element(self
):
597 # Temporary workaround for https://trac.sagemath.org/ticket/28327
598 if self
.is_trivial():
601 s
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
602 return s
.random_element()
607 Return the rank of this EJA.
611 The author knows of no algorithm to compute the rank of an EJA
612 where only the multiplication table is known. In lieu of one, we
613 require the rank to be specified when the algebra is created,
614 and simply pass along that number here.
618 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
619 ....: RealSymmetricEJA,
620 ....: ComplexHermitianEJA,
621 ....: QuaternionHermitianEJA,
626 The rank of the Jordan spin algebra is always two::
628 sage: JordanSpinEJA(2).rank()
630 sage: JordanSpinEJA(3).rank()
632 sage: JordanSpinEJA(4).rank()
635 The rank of the `n`-by-`n` Hermitian real, complex, or
636 quaternion matrices is `n`::
638 sage: RealSymmetricEJA(2).rank()
640 sage: ComplexHermitianEJA(2).rank()
642 sage: QuaternionHermitianEJA(2).rank()
644 sage: RealSymmetricEJA(5).rank()
646 sage: ComplexHermitianEJA(5).rank()
648 sage: QuaternionHermitianEJA(5).rank()
653 Ensure that every EJA that we know how to construct has a
654 positive integer rank::
656 sage: set_random_seed()
657 sage: r = random_eja().rank()
658 sage: r in ZZ and r > 0
665 def vector_space(self
):
667 Return the vector space that underlies this algebra.
671 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
675 sage: J = RealSymmetricEJA(2)
676 sage: J.vector_space()
677 Vector space of dimension 3 over Rational Field
680 return self
.zero().to_vector().parent().ambient_vector_space()
683 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
686 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
688 Return the Euclidean Jordan Algebra corresponding to the set
689 `R^n` under the Hadamard product.
691 Note: this is nothing more than the Cartesian product of ``n``
692 copies of the spin algebra. Once Cartesian product algebras
693 are implemented, this can go.
697 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
701 This multiplication table can be verified by hand::
703 sage: J = RealCartesianProductEJA(3)
704 sage: e0,e1,e2 = J.gens()
720 We can change the generator prefix::
722 sage: RealCartesianProductEJA(3, prefix='r').gens()
725 Our inner product satisfies the Jordan axiom::
727 sage: set_random_seed()
728 sage: n = ZZ.random_element(1,5)
729 sage: J = RealCartesianProductEJA(n)
730 sage: x = J.random_element()
731 sage: y = J.random_element()
732 sage: z = J.random_element()
733 sage: (x*y).inner_product(z) == y.inner_product(x*z)
737 def __init__(self
, n
, field
=QQ
, **kwargs
):
738 V
= VectorSpace(field
, n
)
739 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
742 fdeja
= super(RealCartesianProductEJA
, self
)
743 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
745 def inner_product(self
, x
, y
):
746 return _usual_ip(x
,y
)
751 Return a "random" finite-dimensional Euclidean Jordan Algebra.
755 For now, we choose a random natural number ``n`` (greater than zero)
756 and then give you back one of the following:
758 * The cartesian product of the rational numbers ``n`` times; this is
759 ``QQ^n`` with the Hadamard product.
761 * The Jordan spin algebra on ``QQ^n``.
763 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
766 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
767 in the space of ``2n``-by-``2n`` real symmetric matrices.
769 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
770 in the space of ``4n``-by-``4n`` real symmetric matrices.
772 Later this might be extended to return Cartesian products of the
777 sage: from mjo.eja.eja_algebra import random_eja
782 Euclidean Jordan algebra of dimension...
786 # The max_n component lets us choose different upper bounds on the
787 # value "n" that gets passed to the constructor. This is needed
788 # because e.g. R^{10} is reasonable to test, while the Hermitian
789 # 10-by-10 quaternion matrices are not.
790 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
792 (RealSymmetricEJA
, 5),
793 (ComplexHermitianEJA
, 4),
794 (QuaternionHermitianEJA
, 3)])
795 n
= ZZ
.random_element(1, max_n
)
796 return constructor(n
, field
=QQ
)
800 def _real_symmetric_basis(n
, field
):
802 Return a basis for the space of real symmetric n-by-n matrices.
806 sage: from mjo.eja.eja_algebra import _real_symmetric_basis
810 sage: set_random_seed()
811 sage: n = ZZ.random_element(1,5)
812 sage: B = _real_symmetric_basis(n, QQbar)
813 sage: all( M.is_symmetric() for M in B)
817 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
821 for j
in xrange(i
+1):
822 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
826 # Beware, orthogonal but not normalized!
827 Sij
= Eij
+ Eij
.transpose()
832 def _complex_hermitian_basis(n
, field
):
834 Returns a basis for the space of complex Hermitian n-by-n matrices.
838 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
842 sage: set_random_seed()
843 sage: n = ZZ.random_element(1,5)
844 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
848 F
= QuadraticField(-1, 'I')
851 # This is like the symmetric case, but we need to be careful:
853 # * We want conjugate-symmetry, not just symmetry.
854 # * The diagonal will (as a result) be real.
858 for j
in xrange(i
+1):
859 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
861 Sij
= _embed_complex_matrix(Eij
)
864 # Beware, orthogonal but not normalized! The second one
865 # has a minus because it's conjugated.
866 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
868 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
873 def _quaternion_hermitian_basis(n
, field
):
875 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
879 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
883 sage: set_random_seed()
884 sage: n = ZZ.random_element(1,5)
885 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
889 Q
= QuaternionAlgebra(QQ
,-1,-1)
892 # This is like the symmetric case, but we need to be careful:
894 # * We want conjugate-symmetry, not just symmetry.
895 # * The diagonal will (as a result) be real.
899 for j
in xrange(i
+1):
900 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
902 Sij
= _embed_quaternion_matrix(Eij
)
905 # Beware, orthogonal but not normalized! The second,
906 # third, and fourth ones have a minus because they're
908 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
910 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
912 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
914 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
920 def _multiplication_table_from_matrix_basis(basis
):
922 At least three of the five simple Euclidean Jordan algebras have the
923 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
924 multiplication on the right is matrix multiplication. Given a basis
925 for the underlying matrix space, this function returns a
926 multiplication table (obtained by looping through the basis
927 elements) for an algebra of those matrices.
929 # In S^2, for example, we nominally have four coordinates even
930 # though the space is of dimension three only. The vector space V
931 # is supposed to hold the entire long vector, and the subspace W
932 # of V will be spanned by the vectors that arise from symmetric
933 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
934 field
= basis
[0].base_ring()
935 dimension
= basis
[0].nrows()
937 V
= VectorSpace(field
, dimension
**2)
938 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
940 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
943 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
944 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
949 def _embed_complex_matrix(M
):
951 Embed the n-by-n complex matrix ``M`` into the space of real
952 matrices of size 2n-by-2n via the map the sends each entry `z = a +
953 bi` to the block matrix ``[[a,b],[-b,a]]``.
957 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
961 sage: F = QuadraticField(-1,'i')
962 sage: x1 = F(4 - 2*i)
963 sage: x2 = F(1 + 2*i)
966 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
967 sage: _embed_complex_matrix(M)
976 Embedding is a homomorphism (isomorphism, in fact)::
978 sage: set_random_seed()
979 sage: n = ZZ.random_element(5)
980 sage: F = QuadraticField(-1, 'i')
981 sage: X = random_matrix(F, n)
982 sage: Y = random_matrix(F, n)
983 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
984 sage: expected = _embed_complex_matrix(X*Y)
985 sage: actual == expected
991 raise ValueError("the matrix 'M' must be square")
992 field
= M
.base_ring()
997 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
999 # We can drop the imaginaries here.
1000 return matrix
.block(field
.base_ring(), n
, blocks
)
1003 def _unembed_complex_matrix(M
):
1005 The inverse of _embed_complex_matrix().
1009 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
1010 ....: _unembed_complex_matrix)
1014 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1015 ....: [-2, 1, -4, 3],
1016 ....: [ 9, 10, 11, 12],
1017 ....: [-10, 9, -12, 11] ])
1018 sage: _unembed_complex_matrix(A)
1020 [ 10*i + 9 12*i + 11]
1024 Unembedding is the inverse of embedding::
1026 sage: set_random_seed()
1027 sage: F = QuadraticField(-1, 'i')
1028 sage: M = random_matrix(F, 3)
1029 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1035 raise ValueError("the matrix 'M' must be square")
1036 if not n
.mod(2).is_zero():
1037 raise ValueError("the matrix 'M' must be a complex embedding")
1039 F
= QuadraticField(-1, 'i')
1042 # Go top-left to bottom-right (reading order), converting every
1043 # 2-by-2 block we see to a single complex element.
1045 for k
in xrange(n
/2):
1046 for j
in xrange(n
/2):
1047 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1048 if submat
[0,0] != submat
[1,1]:
1049 raise ValueError('bad on-diagonal submatrix')
1050 if submat
[0,1] != -submat
[1,0]:
1051 raise ValueError('bad off-diagonal submatrix')
1052 z
= submat
[0,0] + submat
[0,1]*i
1055 return matrix(F
, n
/2, elements
)
1058 def _embed_quaternion_matrix(M
):
1060 Embed the n-by-n quaternion matrix ``M`` into the space of real
1061 matrices of size 4n-by-4n by first sending each quaternion entry
1062 `z = a + bi + cj + dk` to the block-complex matrix
1063 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1068 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
1072 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1073 sage: i,j,k = Q.gens()
1074 sage: x = 1 + 2*i + 3*j + 4*k
1075 sage: M = matrix(Q, 1, [[x]])
1076 sage: _embed_quaternion_matrix(M)
1082 Embedding is a homomorphism (isomorphism, in fact)::
1084 sage: set_random_seed()
1085 sage: n = ZZ.random_element(5)
1086 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1087 sage: X = random_matrix(Q, n)
1088 sage: Y = random_matrix(Q, n)
1089 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1090 sage: expected = _embed_quaternion_matrix(X*Y)
1091 sage: actual == expected
1095 quaternions
= M
.base_ring()
1098 raise ValueError("the matrix 'M' must be square")
1100 F
= QuadraticField(-1, 'i')
1105 t
= z
.coefficient_tuple()
1110 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1111 [-c
+ d
*i
, a
- b
*i
]])
1112 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1114 # We should have real entries by now, so use the realest field
1115 # we've got for the return value.
1116 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1119 def _unembed_quaternion_matrix(M
):
1121 The inverse of _embed_quaternion_matrix().
1125 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1126 ....: _unembed_quaternion_matrix)
1130 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1131 ....: [-2, 1, -4, 3],
1132 ....: [-3, 4, 1, -2],
1133 ....: [-4, -3, 2, 1]])
1134 sage: _unembed_quaternion_matrix(M)
1135 [1 + 2*i + 3*j + 4*k]
1139 Unembedding is the inverse of embedding::
1141 sage: set_random_seed()
1142 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1143 sage: M = random_matrix(Q, 3)
1144 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1150 raise ValueError("the matrix 'M' must be square")
1151 if not n
.mod(4).is_zero():
1152 raise ValueError("the matrix 'M' must be a complex embedding")
1154 Q
= QuaternionAlgebra(QQ
,-1,-1)
1157 # Go top-left to bottom-right (reading order), converting every
1158 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1161 for l
in xrange(n
/4):
1162 for m
in xrange(n
/4):
1163 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1164 if submat
[0,0] != submat
[1,1].conjugate():
1165 raise ValueError('bad on-diagonal submatrix')
1166 if submat
[0,1] != -submat
[1,0].conjugate():
1167 raise ValueError('bad off-diagonal submatrix')
1168 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1169 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1172 return matrix(Q
, n
/4, elements
)
1175 # The usual inner product on R^n.
1177 return x
.to_vector().inner_product(y
.to_vector())
1179 # The inner product used for the real symmetric simple EJA.
1180 # We keep it as a separate function because e.g. the complex
1181 # algebra uses the same inner product, except divided by 2.
1182 def _matrix_ip(X
,Y
):
1183 X_mat
= X
.natural_representation()
1184 Y_mat
= Y
.natural_representation()
1185 return (X_mat
*Y_mat
).trace()
1188 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1190 The rank-n simple EJA consisting of real symmetric n-by-n
1191 matrices, the usual symmetric Jordan product, and the trace inner
1192 product. It has dimension `(n^2 + n)/2` over the reals.
1196 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1200 sage: J = RealSymmetricEJA(2)
1201 sage: e0, e1, e2 = J.gens()
1211 The dimension of this algebra is `(n^2 + n) / 2`::
1213 sage: set_random_seed()
1214 sage: n = ZZ.random_element(1,5)
1215 sage: J = RealSymmetricEJA(n)
1216 sage: J.dimension() == (n^2 + n)/2
1219 The Jordan multiplication is what we think it is::
1221 sage: set_random_seed()
1222 sage: n = ZZ.random_element(1,5)
1223 sage: J = RealSymmetricEJA(n)
1224 sage: x = J.random_element()
1225 sage: y = J.random_element()
1226 sage: actual = (x*y).natural_representation()
1227 sage: X = x.natural_representation()
1228 sage: Y = y.natural_representation()
1229 sage: expected = (X*Y + Y*X)/2
1230 sage: actual == expected
1232 sage: J(expected) == x*y
1235 We can change the generator prefix::
1237 sage: RealSymmetricEJA(3, prefix='q').gens()
1238 (q0, q1, q2, q3, q4, q5)
1240 Our inner product satisfies the Jordan axiom::
1242 sage: set_random_seed()
1243 sage: n = ZZ.random_element(1,5)
1244 sage: J = RealSymmetricEJA(n)
1245 sage: x = J.random_element()
1246 sage: y = J.random_element()
1247 sage: z = J.random_element()
1248 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1252 def __init__(self
, n
, field
=QQ
, **kwargs
):
1253 S
= _real_symmetric_basis(n
, field
)
1254 Qs
= _multiplication_table_from_matrix_basis(S
)
1256 fdeja
= super(RealSymmetricEJA
, self
)
1257 return fdeja
.__init
__(field
,
1263 def inner_product(self
, x
, y
):
1264 return _matrix_ip(x
,y
)
1267 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1269 The rank-n simple EJA consisting of complex Hermitian n-by-n
1270 matrices over the real numbers, the usual symmetric Jordan product,
1271 and the real-part-of-trace inner product. It has dimension `n^2` over
1276 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1280 The dimension of this algebra is `n^2`::
1282 sage: set_random_seed()
1283 sage: n = ZZ.random_element(1,5)
1284 sage: J = ComplexHermitianEJA(n)
1285 sage: J.dimension() == n^2
1288 The Jordan multiplication is what we think it is::
1290 sage: set_random_seed()
1291 sage: n = ZZ.random_element(1,5)
1292 sage: J = ComplexHermitianEJA(n)
1293 sage: x = J.random_element()
1294 sage: y = J.random_element()
1295 sage: actual = (x*y).natural_representation()
1296 sage: X = x.natural_representation()
1297 sage: Y = y.natural_representation()
1298 sage: expected = (X*Y + Y*X)/2
1299 sage: actual == expected
1301 sage: J(expected) == x*y
1304 We can change the generator prefix::
1306 sage: ComplexHermitianEJA(2, prefix='z').gens()
1309 Our inner product satisfies the Jordan axiom::
1311 sage: set_random_seed()
1312 sage: n = ZZ.random_element(1,5)
1313 sage: J = ComplexHermitianEJA(n)
1314 sage: x = J.random_element()
1315 sage: y = J.random_element()
1316 sage: z = J.random_element()
1317 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1321 def __init__(self
, n
, field
=QQ
, **kwargs
):
1322 S
= _complex_hermitian_basis(n
, field
)
1323 Qs
= _multiplication_table_from_matrix_basis(S
)
1325 fdeja
= super(ComplexHermitianEJA
, self
)
1326 return fdeja
.__init
__(field
,
1333 def inner_product(self
, x
, y
):
1334 # Since a+bi on the diagonal is represented as
1339 # we'll double-count the "a" entries if we take the trace of
1341 return _matrix_ip(x
,y
)/2
1344 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1346 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1347 matrices, the usual symmetric Jordan product, and the
1348 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1353 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1357 The dimension of this algebra is `n^2`::
1359 sage: set_random_seed()
1360 sage: n = ZZ.random_element(1,5)
1361 sage: J = QuaternionHermitianEJA(n)
1362 sage: J.dimension() == 2*(n^2) - n
1365 The Jordan multiplication is what we think it is::
1367 sage: set_random_seed()
1368 sage: n = ZZ.random_element(1,5)
1369 sage: J = QuaternionHermitianEJA(n)
1370 sage: x = J.random_element()
1371 sage: y = J.random_element()
1372 sage: actual = (x*y).natural_representation()
1373 sage: X = x.natural_representation()
1374 sage: Y = y.natural_representation()
1375 sage: expected = (X*Y + Y*X)/2
1376 sage: actual == expected
1378 sage: J(expected) == x*y
1381 We can change the generator prefix::
1383 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1384 (a0, a1, a2, a3, a4, a5)
1386 Our inner product satisfies the Jordan axiom::
1388 sage: set_random_seed()
1389 sage: n = ZZ.random_element(1,5)
1390 sage: J = QuaternionHermitianEJA(n)
1391 sage: x = J.random_element()
1392 sage: y = J.random_element()
1393 sage: z = J.random_element()
1394 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1398 def __init__(self
, n
, field
=QQ
, **kwargs
):
1399 S
= _quaternion_hermitian_basis(n
, field
)
1400 Qs
= _multiplication_table_from_matrix_basis(S
)
1402 fdeja
= super(QuaternionHermitianEJA
, self
)
1403 return fdeja
.__init
__(field
,
1409 def inner_product(self
, x
, y
):
1410 # Since a+bi+cj+dk on the diagonal is represented as
1412 # a + bi +cj + dk = [ a b c d]
1417 # we'll quadruple-count the "a" entries if we take the trace of
1419 return _matrix_ip(x
,y
)/4
1422 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1424 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1425 with the usual inner product and jordan product ``x*y =
1426 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1431 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1435 This multiplication table can be verified by hand::
1437 sage: J = JordanSpinEJA(4)
1438 sage: e0,e1,e2,e3 = J.gens()
1454 We can change the generator prefix::
1456 sage: JordanSpinEJA(2, prefix='B').gens()
1459 Our inner product satisfies the Jordan axiom::
1461 sage: set_random_seed()
1462 sage: n = ZZ.random_element(1,5)
1463 sage: J = JordanSpinEJA(n)
1464 sage: x = J.random_element()
1465 sage: y = J.random_element()
1466 sage: z = J.random_element()
1467 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1471 def __init__(self
, n
, field
=QQ
, **kwargs
):
1472 V
= VectorSpace(field
, n
)
1473 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1483 z0
= x
.inner_product(y
)
1484 zbar
= y0
*xbar
+ x0
*ybar
1485 z
= V([z0
] + zbar
.list())
1486 mult_table
[i
][j
] = z
1488 # The rank of the spin algebra is two, unless we're in a
1489 # one-dimensional ambient space (because the rank is bounded by
1490 # the ambient dimension).
1491 fdeja
= super(JordanSpinEJA
, self
)
1492 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1494 def inner_product(self
, x
, y
):
1495 return _usual_ip(x
,y
)