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 basis_space
= natural_basis
[0].matrix_space()
139 if elt
not in basis_space
:
140 raise ValueError("not a naturally-represented algebra element")
142 # Thanks for nothing! Matrix spaces aren't vector spaces in
143 # Sage, so we have to figure out its natural-basis coordinates
144 # ourselves. We use the basis space's ring instead of the
145 # element's ring because the basis space might be an algebraic
146 # closure whereas the base ring of the 3-by-3 identity matrix
147 # could be QQ instead of QQbar.
148 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
149 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
150 coords
= W
.coordinate_vector(_mat2vec(elt
))
151 return self
.from_vector(coords
)
156 Return a string representation of ``self``.
160 sage: from mjo.eja.eja_algebra import JordanSpinEJA
164 Ensure that it says what we think it says::
166 sage: JordanSpinEJA(2, field=QQ)
167 Euclidean Jordan algebra of dimension 2 over Rational Field
168 sage: JordanSpinEJA(3, field=RDF)
169 Euclidean Jordan algebra of dimension 3 over Real Double Field
172 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
173 return fmt
.format(self
.dimension(), self
.base_ring())
175 def product_on_basis(self
, i
, j
):
176 return self
._multiplication
_table
[i
][j
]
178 def _a_regular_element(self
):
180 Guess a regular element. Needed to compute the basis for our
181 characteristic polynomial coefficients.
185 sage: from mjo.eja.eja_algebra import random_eja
189 Ensure that this hacky method succeeds for every algebra that we
190 know how to construct::
192 sage: set_random_seed()
193 sage: J = random_eja()
194 sage: J._a_regular_element().is_regular()
199 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
200 if not z
.is_regular():
201 raise ValueError("don't know a regular element")
206 def _charpoly_basis_space(self
):
208 Return the vector space spanned by the basis used in our
209 characteristic polynomial coefficients. This is used not only to
210 compute those coefficients, but also any time we need to
211 evaluate the coefficients (like when we compute the trace or
214 z
= self
._a
_regular
_element
()
215 # Don't use the parent vector space directly here in case this
216 # happens to be a subalgebra. In that case, we would be e.g.
217 # two-dimensional but span_of_basis() would expect three
219 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
220 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
221 V1
= V
.span_of_basis( basis
)
222 b
= (V1
.basis() + V1
.complement().basis())
223 return V
.span_of_basis(b
)
227 def _charpoly_coeff(self
, i
):
229 Return the coefficient polynomial "a_{i}" of this algebra's
230 general characteristic polynomial.
232 Having this be a separate cached method lets us compute and
233 store the trace/determinant (a_{r-1} and a_{0} respectively)
234 separate from the entire characteristic polynomial.
236 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
237 R
= A_of_x
.base_ring()
239 # Guaranteed by theory
242 # Danger: the in-place modification is done for performance
243 # reasons (reconstructing a matrix with huge polynomial
244 # entries is slow), but I don't know how cached_method works,
245 # so it's highly possible that we're modifying some global
246 # list variable by reference, here. In other words, you
247 # probably shouldn't call this method twice on the same
248 # algebra, at the same time, in two threads
249 Ai_orig
= A_of_x
.column(i
)
250 A_of_x
.set_column(i
,xr
)
251 numerator
= A_of_x
.det()
252 A_of_x
.set_column(i
,Ai_orig
)
254 # We're relying on the theory here to ensure that each a_i is
255 # indeed back in R, and the added negative signs are to make
256 # the whole charpoly expression sum to zero.
257 return R(-numerator
/detA
)
261 def _charpoly_matrix_system(self
):
263 Compute the matrix whose entries A_ij are polynomials in
264 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
265 corresponding to `x^r` and the determinent of the matrix A =
266 [A_ij]. In other words, all of the fixed (cachable) data needed
267 to compute the coefficients of the characteristic polynomial.
272 # Turn my vector space into a module so that "vectors" can
273 # have multivatiate polynomial entries.
274 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
275 R
= PolynomialRing(self
.base_ring(), names
)
277 # Using change_ring() on the parent's vector space doesn't work
278 # here because, in a subalgebra, that vector space has a basis
279 # and change_ring() tries to bring the basis along with it. And
280 # that doesn't work unless the new ring is a PID, which it usually
284 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
288 # And figure out the "left multiplication by x" matrix in
291 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
292 for i
in range(n
) ] # don't recompute these!
294 ek
= self
.monomial(k
).to_vector()
296 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
297 for i
in range(n
) ) )
298 Lx
= matrix
.column(R
, lmbx_cols
)
300 # Now we can compute powers of x "symbolically"
301 x_powers
= [self
.one().to_vector(), x
]
302 for d
in range(2, r
+1):
303 x_powers
.append( Lx
*(x_powers
[-1]) )
305 idmat
= matrix
.identity(R
, n
)
307 W
= self
._charpoly
_basis
_space
()
308 W
= W
.change_ring(R
.fraction_field())
310 # Starting with the standard coordinates x = (X1,X2,...,Xn)
311 # and then converting the entries to W-coordinates allows us
312 # to pass in the standard coordinates to the charpoly and get
313 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
316 # W.coordinates(x^2) eval'd at (standard z-coords)
320 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
322 # We want the middle equivalent thing in our matrix, but use
323 # the first equivalent thing instead so that we can pass in
324 # standard coordinates.
325 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
326 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
327 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
328 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
332 def characteristic_polynomial(self
):
334 Return a characteristic polynomial that works for all elements
337 The resulting polynomial has `n+1` variables, where `n` is the
338 dimension of this algebra. The first `n` variables correspond to
339 the coordinates of an algebra element: when evaluated at the
340 coordinates of an algebra element with respect to a certain
341 basis, the result is a univariate polynomial (in the one
342 remaining variable ``t``), namely the characteristic polynomial
347 sage: from mjo.eja.eja_algebra import JordanSpinEJA
351 The characteristic polynomial in the spin algebra is given in
352 Alizadeh, Example 11.11::
354 sage: J = JordanSpinEJA(3)
355 sage: p = J.characteristic_polynomial(); p
356 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
357 sage: xvec = J.one().to_vector()
365 # The list of coefficient polynomials a_1, a_2, ..., a_n.
366 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
368 # We go to a bit of trouble here to reorder the
369 # indeterminates, so that it's easier to evaluate the
370 # characteristic polynomial at x's coordinates and get back
371 # something in terms of t, which is what we want.
373 S
= PolynomialRing(self
.base_ring(),'t')
375 S
= PolynomialRing(S
, R
.variable_names())
378 # Note: all entries past the rth should be zero. The
379 # coefficient of the highest power (x^r) is 1, but it doesn't
380 # appear in the solution vector which contains coefficients
381 # for the other powers (to make them sum to x^r).
383 a
[r
] = 1 # corresponds to x^r
385 # When the rank is equal to the dimension, trying to
386 # assign a[r] goes out-of-bounds.
387 a
.append(1) # corresponds to x^r
389 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
392 def inner_product(self
, x
, y
):
394 The inner product associated with this Euclidean Jordan algebra.
396 Defaults to the trace inner product, but can be overridden by
397 subclasses if they are sure that the necessary properties are
402 sage: from mjo.eja.eja_algebra import random_eja
406 The inner product must satisfy its axiom for this algebra to truly
407 be a Euclidean Jordan Algebra::
409 sage: set_random_seed()
410 sage: J = random_eja()
411 sage: x = J.random_element()
412 sage: y = J.random_element()
413 sage: z = J.random_element()
414 sage: (x*y).inner_product(z) == y.inner_product(x*z)
418 if (not x
in self
) or (not y
in self
):
419 raise TypeError("arguments must live in this algebra")
420 return x
.trace_inner_product(y
)
423 def is_trivial(self
):
425 Return whether or not this algebra is trivial.
427 A trivial algebra contains only the zero element.
431 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
435 sage: J = ComplexHermitianEJA(3)
438 sage: A = J.zero().subalgebra_generated_by()
443 return self
.dimension() == 0
446 def multiplication_table(self
):
448 Return a visual representation of this algebra's multiplication
449 table (on basis elements).
453 sage: from mjo.eja.eja_algebra import JordanSpinEJA
457 sage: J = JordanSpinEJA(4)
458 sage: J.multiplication_table()
459 +----++----+----+----+----+
460 | * || e0 | e1 | e2 | e3 |
461 +====++====+====+====+====+
462 | e0 || e0 | e1 | e2 | e3 |
463 +----++----+----+----+----+
464 | e1 || e1 | e0 | 0 | 0 |
465 +----++----+----+----+----+
466 | e2 || e2 | 0 | e0 | 0 |
467 +----++----+----+----+----+
468 | e3 || e3 | 0 | 0 | e0 |
469 +----++----+----+----+----+
472 M
= list(self
._multiplication
_table
) # copy
473 for i
in range(len(M
)):
474 # M had better be "square"
475 M
[i
] = [self
.monomial(i
)] + M
[i
]
476 M
= [["*"] + list(self
.gens())] + M
477 return table(M
, header_row
=True, header_column
=True, frame
=True)
480 def natural_basis(self
):
482 Return a more-natural representation of this algebra's basis.
484 Every finite-dimensional Euclidean Jordan Algebra is a direct
485 sum of five simple algebras, four of which comprise Hermitian
486 matrices. This method returns the original "natural" basis
487 for our underlying vector space. (Typically, the natural basis
488 is used to construct the multiplication table in the first place.)
490 Note that this will always return a matrix. The standard basis
491 in `R^n` will be returned as `n`-by-`1` column matrices.
495 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
496 ....: RealSymmetricEJA)
500 sage: J = RealSymmetricEJA(2)
502 Finite family {0: e0, 1: e1, 2: e2}
503 sage: J.natural_basis()
505 [1 0] [ 0 1/2*sqrt2] [0 0]
506 [0 0], [1/2*sqrt2 0], [0 1]
511 sage: J = JordanSpinEJA(2)
513 Finite family {0: e0, 1: e1}
514 sage: J.natural_basis()
521 if self
._natural
_basis
is None:
522 M
= self
.natural_basis_space()
523 return tuple( M(b
.to_vector()) for b
in self
.basis() )
525 return self
._natural
_basis
528 def natural_basis_space(self
):
530 Return the matrix space in which this algebra's natural basis
533 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
534 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
536 return self
._natural
_basis
[0].matrix_space()
542 Return the unit element of this algebra.
546 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
551 sage: J = RealCartesianProductEJA(5)
553 e0 + e1 + e2 + e3 + e4
557 The identity element acts like the identity::
559 sage: set_random_seed()
560 sage: J = random_eja()
561 sage: x = J.random_element()
562 sage: J.one()*x == x and x*J.one() == x
565 The matrix of the unit element's operator is the identity::
567 sage: set_random_seed()
568 sage: J = random_eja()
569 sage: actual = J.one().operator().matrix()
570 sage: expected = matrix.identity(J.base_ring(), J.dimension())
571 sage: actual == expected
575 # We can brute-force compute the matrices of the operators
576 # that correspond to the basis elements of this algebra.
577 # If some linear combination of those basis elements is the
578 # algebra identity, then the same linear combination of
579 # their matrices has to be the identity matrix.
581 # Of course, matrices aren't vectors in sage, so we have to
582 # appeal to the "long vectors" isometry.
583 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
585 # Now we use basis linear algebra to find the coefficients,
586 # of the matrices-as-vectors-linear-combination, which should
587 # work for the original algebra basis too.
588 A
= matrix
.column(self
.base_ring(), oper_vecs
)
590 # We used the isometry on the left-hand side already, but we
591 # still need to do it for the right-hand side. Recall that we
592 # wanted something that summed to the identity matrix.
593 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
595 # Now if there's an identity element in the algebra, this should work.
596 coeffs
= A
.solve_right(b
)
597 return self
.linear_combination(zip(self
.gens(), coeffs
))
600 def random_element(self
):
601 # Temporary workaround for https://trac.sagemath.org/ticket/28327
602 if self
.is_trivial():
605 s
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
606 return s
.random_element()
611 Return the rank of this EJA.
615 The author knows of no algorithm to compute the rank of an EJA
616 where only the multiplication table is known. In lieu of one, we
617 require the rank to be specified when the algebra is created,
618 and simply pass along that number here.
622 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
623 ....: RealSymmetricEJA,
624 ....: ComplexHermitianEJA,
625 ....: QuaternionHermitianEJA,
630 The rank of the Jordan spin algebra is always two::
632 sage: JordanSpinEJA(2).rank()
634 sage: JordanSpinEJA(3).rank()
636 sage: JordanSpinEJA(4).rank()
639 The rank of the `n`-by-`n` Hermitian real, complex, or
640 quaternion matrices is `n`::
642 sage: RealSymmetricEJA(2).rank()
644 sage: ComplexHermitianEJA(2).rank()
646 sage: QuaternionHermitianEJA(2).rank()
648 sage: RealSymmetricEJA(5).rank()
650 sage: ComplexHermitianEJA(5).rank()
652 sage: QuaternionHermitianEJA(5).rank()
657 Ensure that every EJA that we know how to construct has a
658 positive integer rank::
660 sage: set_random_seed()
661 sage: r = random_eja().rank()
662 sage: r in ZZ and r > 0
669 def vector_space(self
):
671 Return the vector space that underlies this algebra.
675 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
679 sage: J = RealSymmetricEJA(2)
680 sage: J.vector_space()
681 Vector space of dimension 3 over...
684 return self
.zero().to_vector().parent().ambient_vector_space()
687 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
690 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
692 Return the Euclidean Jordan Algebra corresponding to the set
693 `R^n` under the Hadamard product.
695 Note: this is nothing more than the Cartesian product of ``n``
696 copies of the spin algebra. Once Cartesian product algebras
697 are implemented, this can go.
701 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
705 This multiplication table can be verified by hand::
707 sage: J = RealCartesianProductEJA(3)
708 sage: e0,e1,e2 = J.gens()
724 We can change the generator prefix::
726 sage: RealCartesianProductEJA(3, prefix='r').gens()
729 Our inner product satisfies the Jordan axiom::
731 sage: set_random_seed()
732 sage: n = ZZ.random_element(1,5)
733 sage: J = RealCartesianProductEJA(n)
734 sage: x = J.random_element()
735 sage: y = J.random_element()
736 sage: z = J.random_element()
737 sage: (x*y).inner_product(z) == y.inner_product(x*z)
741 def __init__(self
, n
, field
=QQ
, **kwargs
):
742 V
= VectorSpace(field
, n
)
743 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
746 fdeja
= super(RealCartesianProductEJA
, self
)
747 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
749 def inner_product(self
, x
, y
):
750 return _usual_ip(x
,y
)
755 Return a "random" finite-dimensional Euclidean Jordan Algebra.
759 For now, we choose a random natural number ``n`` (greater than zero)
760 and then give you back one of the following:
762 * The cartesian product of the rational numbers ``n`` times; this is
763 ``QQ^n`` with the Hadamard product.
765 * The Jordan spin algebra on ``QQ^n``.
767 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
770 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
771 in the space of ``2n``-by-``2n`` real symmetric matrices.
773 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
774 in the space of ``4n``-by-``4n`` real symmetric matrices.
776 Later this might be extended to return Cartesian products of the
781 sage: from mjo.eja.eja_algebra import random_eja
786 Euclidean Jordan algebra of dimension...
790 # The max_n component lets us choose different upper bounds on the
791 # value "n" that gets passed to the constructor. This is needed
792 # because e.g. R^{10} is reasonable to test, while the Hermitian
793 # 10-by-10 quaternion matrices are not.
794 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
796 (RealSymmetricEJA
, 5),
797 (ComplexHermitianEJA
, 4),
798 (QuaternionHermitianEJA
, 3)])
799 n
= ZZ
.random_element(1, max_n
)
800 return constructor(n
, field
=QQ
)
804 def _real_symmetric_basis(n
, field
):
806 Return a basis for the space of real symmetric n-by-n matrices.
810 sage: from mjo.eja.eja_algebra import _real_symmetric_basis
814 sage: set_random_seed()
815 sage: n = ZZ.random_element(1,5)
816 sage: B = _real_symmetric_basis(n, QQbar)
817 sage: all( M.is_symmetric() for M in B)
821 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
825 for j
in xrange(i
+1):
826 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
830 Sij
= Eij
+ Eij
.transpose()
832 Sij
= Sij
/ _real_symmetric_matrix_ip(Sij
,Sij
).sqrt()
837 def _complex_hermitian_basis(n
, field
):
839 Returns a basis for the space of complex Hermitian n-by-n matrices.
843 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
847 sage: set_random_seed()
848 sage: n = ZZ.random_element(1,5)
849 sage: B = _complex_hermitian_basis(n, QQ)
850 sage: all( M.is_symmetric() for M in B)
854 F
= QuadraticField(-1, 'I')
857 # This is like the symmetric case, but we need to be careful:
859 # * We want conjugate-symmetry, not just symmetry.
860 # * The diagonal will (as a result) be real.
864 for j
in xrange(i
+1):
865 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
867 Sij
= _embed_complex_matrix(Eij
)
870 # Beware, orthogonal but not normalized! The second one
871 # has a minus because it's conjugated.
872 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
874 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
879 def _quaternion_hermitian_basis(n
, field
):
881 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
885 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
889 sage: set_random_seed()
890 sage: n = ZZ.random_element(1,5)
891 sage: B = _quaternion_hermitian_basis(n, QQbar)
892 sage: all( M.is_symmetric() for M in B )
896 Q
= QuaternionAlgebra(QQ
,-1,-1)
899 # This is like the symmetric case, but we need to be careful:
901 # * We want conjugate-symmetry, not just symmetry.
902 # * The diagonal will (as a result) be real.
906 for j
in xrange(i
+1):
907 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
909 Sij
= _embed_quaternion_matrix(Eij
)
912 # Beware, orthogonal but not normalized! The second,
913 # third, and fourth ones have a minus because they're
915 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
917 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
919 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
921 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
927 def _multiplication_table_from_matrix_basis(basis
):
929 At least three of the five simple Euclidean Jordan algebras have the
930 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
931 multiplication on the right is matrix multiplication. Given a basis
932 for the underlying matrix space, this function returns a
933 multiplication table (obtained by looping through the basis
934 elements) for an algebra of those matrices.
936 # In S^2, for example, we nominally have four coordinates even
937 # though the space is of dimension three only. The vector space V
938 # is supposed to hold the entire long vector, and the subspace W
939 # of V will be spanned by the vectors that arise from symmetric
940 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
941 field
= basis
[0].base_ring()
942 dimension
= basis
[0].nrows()
944 V
= VectorSpace(field
, dimension
**2)
945 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
947 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
950 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
951 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
956 def _embed_complex_matrix(M
):
958 Embed the n-by-n complex matrix ``M`` into the space of real
959 matrices of size 2n-by-2n via the map the sends each entry `z = a +
960 bi` to the block matrix ``[[a,b],[-b,a]]``.
964 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
968 sage: F = QuadraticField(-1,'i')
969 sage: x1 = F(4 - 2*i)
970 sage: x2 = F(1 + 2*i)
973 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
974 sage: _embed_complex_matrix(M)
983 Embedding is a homomorphism (isomorphism, in fact)::
985 sage: set_random_seed()
986 sage: n = ZZ.random_element(5)
987 sage: F = QuadraticField(-1, 'i')
988 sage: X = random_matrix(F, n)
989 sage: Y = random_matrix(F, n)
990 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
991 sage: expected = _embed_complex_matrix(X*Y)
992 sage: actual == expected
998 raise ValueError("the matrix 'M' must be square")
999 field
= M
.base_ring()
1004 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1006 # We can drop the imaginaries here.
1007 return matrix
.block(field
.base_ring(), n
, blocks
)
1010 def _unembed_complex_matrix(M
):
1012 The inverse of _embed_complex_matrix().
1016 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
1017 ....: _unembed_complex_matrix)
1021 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1022 ....: [-2, 1, -4, 3],
1023 ....: [ 9, 10, 11, 12],
1024 ....: [-10, 9, -12, 11] ])
1025 sage: _unembed_complex_matrix(A)
1027 [ 10*i + 9 12*i + 11]
1031 Unembedding is the inverse of embedding::
1033 sage: set_random_seed()
1034 sage: F = QuadraticField(-1, 'i')
1035 sage: M = random_matrix(F, 3)
1036 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1042 raise ValueError("the matrix 'M' must be square")
1043 if not n
.mod(2).is_zero():
1044 raise ValueError("the matrix 'M' must be a complex embedding")
1046 F
= QuadraticField(-1, 'i')
1049 # Go top-left to bottom-right (reading order), converting every
1050 # 2-by-2 block we see to a single complex element.
1052 for k
in xrange(n
/2):
1053 for j
in xrange(n
/2):
1054 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1055 if submat
[0,0] != submat
[1,1]:
1056 raise ValueError('bad on-diagonal submatrix')
1057 if submat
[0,1] != -submat
[1,0]:
1058 raise ValueError('bad off-diagonal submatrix')
1059 z
= submat
[0,0] + submat
[0,1]*i
1062 return matrix(F
, n
/2, elements
)
1065 def _embed_quaternion_matrix(M
):
1067 Embed the n-by-n quaternion matrix ``M`` into the space of real
1068 matrices of size 4n-by-4n by first sending each quaternion entry
1069 `z = a + bi + cj + dk` to the block-complex matrix
1070 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1075 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
1079 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1080 sage: i,j,k = Q.gens()
1081 sage: x = 1 + 2*i + 3*j + 4*k
1082 sage: M = matrix(Q, 1, [[x]])
1083 sage: _embed_quaternion_matrix(M)
1089 Embedding is a homomorphism (isomorphism, in fact)::
1091 sage: set_random_seed()
1092 sage: n = ZZ.random_element(5)
1093 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1094 sage: X = random_matrix(Q, n)
1095 sage: Y = random_matrix(Q, n)
1096 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1097 sage: expected = _embed_quaternion_matrix(X*Y)
1098 sage: actual == expected
1102 quaternions
= M
.base_ring()
1105 raise ValueError("the matrix 'M' must be square")
1107 F
= QuadraticField(-1, 'i')
1112 t
= z
.coefficient_tuple()
1117 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1118 [-c
+ d
*i
, a
- b
*i
]])
1119 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1121 # We should have real entries by now, so use the realest field
1122 # we've got for the return value.
1123 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1126 def _unembed_quaternion_matrix(M
):
1128 The inverse of _embed_quaternion_matrix().
1132 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1133 ....: _unembed_quaternion_matrix)
1137 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1138 ....: [-2, 1, -4, 3],
1139 ....: [-3, 4, 1, -2],
1140 ....: [-4, -3, 2, 1]])
1141 sage: _unembed_quaternion_matrix(M)
1142 [1 + 2*i + 3*j + 4*k]
1146 Unembedding is the inverse of embedding::
1148 sage: set_random_seed()
1149 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1150 sage: M = random_matrix(Q, 3)
1151 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1157 raise ValueError("the matrix 'M' must be square")
1158 if not n
.mod(4).is_zero():
1159 raise ValueError("the matrix 'M' must be a complex embedding")
1161 Q
= QuaternionAlgebra(QQ
,-1,-1)
1164 # Go top-left to bottom-right (reading order), converting every
1165 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1168 for l
in xrange(n
/4):
1169 for m
in xrange(n
/4):
1170 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1171 if submat
[0,0] != submat
[1,1].conjugate():
1172 raise ValueError('bad on-diagonal submatrix')
1173 if submat
[0,1] != -submat
[1,0].conjugate():
1174 raise ValueError('bad off-diagonal submatrix')
1175 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1176 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1179 return matrix(Q
, n
/4, elements
)
1182 # The usual inner product on R^n.
1184 return x
.to_vector().inner_product(y
.to_vector())
1186 # The inner product used for the real symmetric simple EJA.
1187 # We keep it as a separate function because e.g. the complex
1188 # algebra uses the same inner product, except divided by 2.
1189 def _matrix_ip(X
,Y
):
1190 X_mat
= X
.natural_representation()
1191 Y_mat
= Y
.natural_representation()
1192 return (X_mat
*Y_mat
).trace()
1194 def _real_symmetric_matrix_ip(X
,Y
):
1195 return (X
*Y
).trace()
1198 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1200 The rank-n simple EJA consisting of real symmetric n-by-n
1201 matrices, the usual symmetric Jordan product, and the trace inner
1202 product. It has dimension `(n^2 + n)/2` over the reals.
1206 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1210 sage: J = RealSymmetricEJA(2)
1211 sage: e0, e1, e2 = J.gens()
1221 The dimension of this algebra is `(n^2 + n) / 2`::
1223 sage: set_random_seed()
1224 sage: n = ZZ.random_element(1,5)
1225 sage: J = RealSymmetricEJA(n)
1226 sage: J.dimension() == (n^2 + n)/2
1229 The Jordan multiplication is what we think it is::
1231 sage: set_random_seed()
1232 sage: n = ZZ.random_element(1,5)
1233 sage: J = RealSymmetricEJA(n)
1234 sage: x = J.random_element()
1235 sage: y = J.random_element()
1236 sage: actual = (x*y).natural_representation()
1237 sage: X = x.natural_representation()
1238 sage: Y = y.natural_representation()
1239 sage: expected = (X*Y + Y*X)/2
1240 sage: actual == expected
1242 sage: J(expected) == x*y
1245 We can change the generator prefix::
1247 sage: RealSymmetricEJA(3, prefix='q').gens()
1248 (q0, q1, q2, q3, q4, q5)
1250 Our inner product satisfies the Jordan axiom::
1252 sage: set_random_seed()
1253 sage: n = ZZ.random_element(1,5)
1254 sage: J = RealSymmetricEJA(n)
1255 sage: x = J.random_element()
1256 sage: y = J.random_element()
1257 sage: z = J.random_element()
1258 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1261 Our basis is normalized with respect to the natural inner product::
1263 sage: set_random_seed()
1264 sage: n = ZZ.random_element(1,5)
1265 sage: J = RealSymmetricEJA(n)
1266 sage: all( b.norm() == 1 for b in J.gens() )
1269 Left-multiplication operators are symmetric because they satisfy
1272 sage: set_random_seed()
1273 sage: n = ZZ.random_element(1,5)
1274 sage: x = RealSymmetricEJA(n).random_element()
1275 sage: x.operator().matrix().is_symmetric()
1279 def __init__(self
, n
, field
=QQ
, **kwargs
):
1280 if n
> 1 and field
is QQ
:
1281 # We'll need sqrt(2) to normalize the basis, and this
1282 # winds up in the multiplication table, so the whole
1283 # algebra needs to be over the field extension.
1284 field
= QuadraticField(2, 'sqrt2')
1286 S
= _real_symmetric_basis(n
, field
)
1287 Qs
= _multiplication_table_from_matrix_basis(S
)
1289 fdeja
= super(RealSymmetricEJA
, self
)
1290 return fdeja
.__init
__(field
,
1296 def inner_product(self
, x
, y
):
1297 X
= x
.natural_representation()
1298 Y
= y
.natural_representation()
1299 return _real_symmetric_matrix_ip(X
,Y
)
1302 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1304 The rank-n simple EJA consisting of complex Hermitian n-by-n
1305 matrices over the real numbers, the usual symmetric Jordan product,
1306 and the real-part-of-trace inner product. It has dimension `n^2` over
1311 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1315 The dimension of this algebra is `n^2`::
1317 sage: set_random_seed()
1318 sage: n = ZZ.random_element(1,5)
1319 sage: J = ComplexHermitianEJA(n)
1320 sage: J.dimension() == n^2
1323 The Jordan multiplication is what we think it is::
1325 sage: set_random_seed()
1326 sage: n = ZZ.random_element(1,5)
1327 sage: J = ComplexHermitianEJA(n)
1328 sage: x = J.random_element()
1329 sage: y = J.random_element()
1330 sage: actual = (x*y).natural_representation()
1331 sage: X = x.natural_representation()
1332 sage: Y = y.natural_representation()
1333 sage: expected = (X*Y + Y*X)/2
1334 sage: actual == expected
1336 sage: J(expected) == x*y
1339 We can change the generator prefix::
1341 sage: ComplexHermitianEJA(2, prefix='z').gens()
1344 Our inner product satisfies the Jordan axiom::
1346 sage: set_random_seed()
1347 sage: n = ZZ.random_element(1,5)
1348 sage: J = ComplexHermitianEJA(n)
1349 sage: x = J.random_element()
1350 sage: y = J.random_element()
1351 sage: z = J.random_element()
1352 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1356 def __init__(self
, n
, field
=QQ
, **kwargs
):
1357 S
= _complex_hermitian_basis(n
, field
)
1358 Qs
= _multiplication_table_from_matrix_basis(S
)
1360 fdeja
= super(ComplexHermitianEJA
, self
)
1361 return fdeja
.__init
__(field
,
1368 def inner_product(self
, x
, y
):
1369 # Since a+bi on the diagonal is represented as
1374 # we'll double-count the "a" entries if we take the trace of
1376 return _matrix_ip(x
,y
)/2
1379 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1381 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1382 matrices, the usual symmetric Jordan product, and the
1383 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1388 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1392 The dimension of this algebra is `n^2`::
1394 sage: set_random_seed()
1395 sage: n = ZZ.random_element(1,5)
1396 sage: J = QuaternionHermitianEJA(n)
1397 sage: J.dimension() == 2*(n^2) - n
1400 The Jordan multiplication is what we think it is::
1402 sage: set_random_seed()
1403 sage: n = ZZ.random_element(1,5)
1404 sage: J = QuaternionHermitianEJA(n)
1405 sage: x = J.random_element()
1406 sage: y = J.random_element()
1407 sage: actual = (x*y).natural_representation()
1408 sage: X = x.natural_representation()
1409 sage: Y = y.natural_representation()
1410 sage: expected = (X*Y + Y*X)/2
1411 sage: actual == expected
1413 sage: J(expected) == x*y
1416 We can change the generator prefix::
1418 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1419 (a0, a1, a2, a3, a4, a5)
1421 Our inner product satisfies the Jordan axiom::
1423 sage: set_random_seed()
1424 sage: n = ZZ.random_element(1,5)
1425 sage: J = QuaternionHermitianEJA(n)
1426 sage: x = J.random_element()
1427 sage: y = J.random_element()
1428 sage: z = J.random_element()
1429 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1433 def __init__(self
, n
, field
=QQ
, **kwargs
):
1434 S
= _quaternion_hermitian_basis(n
, field
)
1435 Qs
= _multiplication_table_from_matrix_basis(S
)
1437 fdeja
= super(QuaternionHermitianEJA
, self
)
1438 return fdeja
.__init
__(field
,
1444 def inner_product(self
, x
, y
):
1445 # Since a+bi+cj+dk on the diagonal is represented as
1447 # a + bi +cj + dk = [ a b c d]
1452 # we'll quadruple-count the "a" entries if we take the trace of
1454 return _matrix_ip(x
,y
)/4
1457 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1459 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1460 with the usual inner product and jordan product ``x*y =
1461 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1466 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1470 This multiplication table can be verified by hand::
1472 sage: J = JordanSpinEJA(4)
1473 sage: e0,e1,e2,e3 = J.gens()
1489 We can change the generator prefix::
1491 sage: JordanSpinEJA(2, prefix='B').gens()
1494 Our inner product satisfies the Jordan axiom::
1496 sage: set_random_seed()
1497 sage: n = ZZ.random_element(1,5)
1498 sage: J = JordanSpinEJA(n)
1499 sage: x = J.random_element()
1500 sage: y = J.random_element()
1501 sage: z = J.random_element()
1502 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1506 def __init__(self
, n
, field
=QQ
, **kwargs
):
1507 V
= VectorSpace(field
, n
)
1508 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1518 z0
= x
.inner_product(y
)
1519 zbar
= y0
*xbar
+ x0
*ybar
1520 z
= V([z0
] + zbar
.list())
1521 mult_table
[i
][j
] = z
1523 # The rank of the spin algebra is two, unless we're in a
1524 # one-dimensional ambient space (because the rank is bounded by
1525 # the ambient dimension).
1526 fdeja
= super(JordanSpinEJA
, self
)
1527 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1529 def inner_product(self
, x
, y
):
1530 return _usual_ip(x
,y
)