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 is_trivial(self
):
409 Return whether or not this algebra is trivial.
411 A trivial algebra contains only the zero element.
415 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
419 sage: J = ComplexHermitianEJA(3)
422 sage: A = J.zero().subalgebra_generated_by()
427 return self
.dimension() == 0
430 def multiplication_table(self
):
432 Return a visual representation of this algebra's multiplication
433 table (on basis elements).
437 sage: from mjo.eja.eja_algebra import JordanSpinEJA
441 sage: J = JordanSpinEJA(4)
442 sage: J.multiplication_table()
443 +----++----+----+----+----+
444 | * || e0 | e1 | e2 | e3 |
445 +====++====+====+====+====+
446 | e0 || e0 | e1 | e2 | e3 |
447 +----++----+----+----+----+
448 | e1 || e1 | e0 | 0 | 0 |
449 +----++----+----+----+----+
450 | e2 || e2 | 0 | e0 | 0 |
451 +----++----+----+----+----+
452 | e3 || e3 | 0 | 0 | e0 |
453 +----++----+----+----+----+
456 M
= list(self
._multiplication
_table
) # copy
457 for i
in range(len(M
)):
458 # M had better be "square"
459 M
[i
] = [self
.monomial(i
)] + M
[i
]
460 M
= [["*"] + list(self
.gens())] + M
461 return table(M
, header_row
=True, header_column
=True, frame
=True)
464 def natural_basis(self
):
466 Return a more-natural representation of this algebra's basis.
468 Every finite-dimensional Euclidean Jordan Algebra is a direct
469 sum of five simple algebras, four of which comprise Hermitian
470 matrices. This method returns the original "natural" basis
471 for our underlying vector space. (Typically, the natural basis
472 is used to construct the multiplication table in the first place.)
474 Note that this will always return a matrix. The standard basis
475 in `R^n` will be returned as `n`-by-`1` column matrices.
479 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
480 ....: RealSymmetricEJA)
484 sage: J = RealSymmetricEJA(2)
486 Finite family {0: e0, 1: e1, 2: e2}
487 sage: J.natural_basis()
495 sage: J = JordanSpinEJA(2)
497 Finite family {0: e0, 1: e1}
498 sage: J.natural_basis()
505 if self
._natural
_basis
is None:
506 return tuple( b
.to_vector().column() for b
in self
.basis() )
508 return self
._natural
_basis
514 Return the unit element of this algebra.
518 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
523 sage: J = RealCartesianProductEJA(5)
525 e0 + e1 + e2 + e3 + e4
529 The identity element acts like the identity::
531 sage: set_random_seed()
532 sage: J = random_eja()
533 sage: x = J.random_element()
534 sage: J.one()*x == x and x*J.one() == x
537 The matrix of the unit element's operator is the identity::
539 sage: set_random_seed()
540 sage: J = random_eja()
541 sage: actual = J.one().operator().matrix()
542 sage: expected = matrix.identity(J.base_ring(), J.dimension())
543 sage: actual == expected
547 # We can brute-force compute the matrices of the operators
548 # that correspond to the basis elements of this algebra.
549 # If some linear combination of those basis elements is the
550 # algebra identity, then the same linear combination of
551 # their matrices has to be the identity matrix.
553 # Of course, matrices aren't vectors in sage, so we have to
554 # appeal to the "long vectors" isometry.
555 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
557 # Now we use basis linear algebra to find the coefficients,
558 # of the matrices-as-vectors-linear-combination, which should
559 # work for the original algebra basis too.
560 A
= matrix
.column(self
.base_ring(), oper_vecs
)
562 # We used the isometry on the left-hand side already, but we
563 # still need to do it for the right-hand side. Recall that we
564 # wanted something that summed to the identity matrix.
565 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
567 # Now if there's an identity element in the algebra, this should work.
568 coeffs
= A
.solve_right(b
)
569 return self
.linear_combination(zip(self
.gens(), coeffs
))
572 def random_element(self
):
573 # Temporary workaround for https://trac.sagemath.org/ticket/28327
574 if self
.is_trivial():
577 s
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
578 return s
.random_element()
583 Return the rank of this EJA.
587 The author knows of no algorithm to compute the rank of an EJA
588 where only the multiplication table is known. In lieu of one, we
589 require the rank to be specified when the algebra is created,
590 and simply pass along that number here.
594 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
595 ....: RealSymmetricEJA,
596 ....: ComplexHermitianEJA,
597 ....: QuaternionHermitianEJA,
602 The rank of the Jordan spin algebra is always two::
604 sage: JordanSpinEJA(2).rank()
606 sage: JordanSpinEJA(3).rank()
608 sage: JordanSpinEJA(4).rank()
611 The rank of the `n`-by-`n` Hermitian real, complex, or
612 quaternion matrices is `n`::
614 sage: RealSymmetricEJA(2).rank()
616 sage: ComplexHermitianEJA(2).rank()
618 sage: QuaternionHermitianEJA(2).rank()
620 sage: RealSymmetricEJA(5).rank()
622 sage: ComplexHermitianEJA(5).rank()
624 sage: QuaternionHermitianEJA(5).rank()
629 Ensure that every EJA that we know how to construct has a
630 positive integer rank::
632 sage: set_random_seed()
633 sage: r = random_eja().rank()
634 sage: r in ZZ and r > 0
641 def vector_space(self
):
643 Return the vector space that underlies this algebra.
647 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
651 sage: J = RealSymmetricEJA(2)
652 sage: J.vector_space()
653 Vector space of dimension 3 over Rational Field
656 return self
.zero().to_vector().parent().ambient_vector_space()
659 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
662 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
664 Return the Euclidean Jordan Algebra corresponding to the set
665 `R^n` under the Hadamard product.
667 Note: this is nothing more than the Cartesian product of ``n``
668 copies of the spin algebra. Once Cartesian product algebras
669 are implemented, this can go.
673 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
677 This multiplication table can be verified by hand::
679 sage: J = RealCartesianProductEJA(3)
680 sage: e0,e1,e2 = J.gens()
696 We can change the generator prefix::
698 sage: RealCartesianProductEJA(3, prefix='r').gens()
702 def __init__(self
, n
, field
=QQ
, **kwargs
):
703 V
= VectorSpace(field
, n
)
704 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
707 fdeja
= super(RealCartesianProductEJA
, self
)
708 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
710 def inner_product(self
, x
, y
):
711 return _usual_ip(x
,y
)
716 Return a "random" finite-dimensional Euclidean Jordan Algebra.
720 For now, we choose a random natural number ``n`` (greater than zero)
721 and then give you back one of the following:
723 * The cartesian product of the rational numbers ``n`` times; this is
724 ``QQ^n`` with the Hadamard product.
726 * The Jordan spin algebra on ``QQ^n``.
728 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
731 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
732 in the space of ``2n``-by-``2n`` real symmetric matrices.
734 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
735 in the space of ``4n``-by-``4n`` real symmetric matrices.
737 Later this might be extended to return Cartesian products of the
742 sage: from mjo.eja.eja_algebra import random_eja
747 Euclidean Jordan algebra of dimension...
751 # The max_n component lets us choose different upper bounds on the
752 # value "n" that gets passed to the constructor. This is needed
753 # because e.g. R^{10} is reasonable to test, while the Hermitian
754 # 10-by-10 quaternion matrices are not.
755 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
757 (RealSymmetricEJA
, 5),
758 (ComplexHermitianEJA
, 4),
759 (QuaternionHermitianEJA
, 3)])
760 n
= ZZ
.random_element(1, max_n
)
761 return constructor(n
, field
=QQ
)
765 def _real_symmetric_basis(n
, field
=QQ
):
767 Return a basis for the space of real symmetric n-by-n matrices.
769 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
773 for j
in xrange(i
+1):
774 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
778 # Beware, orthogonal but not normalized!
779 Sij
= Eij
+ Eij
.transpose()
784 def _complex_hermitian_basis(n
, field
=QQ
):
786 Returns a basis for the space of complex Hermitian n-by-n matrices.
790 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
794 sage: set_random_seed()
795 sage: n = ZZ.random_element(1,5)
796 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
800 F
= QuadraticField(-1, 'I')
803 # This is like the symmetric case, but we need to be careful:
805 # * We want conjugate-symmetry, not just symmetry.
806 # * The diagonal will (as a result) be real.
810 for j
in xrange(i
+1):
811 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
813 Sij
= _embed_complex_matrix(Eij
)
816 # Beware, orthogonal but not normalized! The second one
817 # has a minus because it's conjugated.
818 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
820 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
825 def _quaternion_hermitian_basis(n
, field
=QQ
):
827 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
831 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
835 sage: set_random_seed()
836 sage: n = ZZ.random_element(1,5)
837 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
841 Q
= QuaternionAlgebra(QQ
,-1,-1)
844 # This is like the symmetric case, but we need to be careful:
846 # * We want conjugate-symmetry, not just symmetry.
847 # * The diagonal will (as a result) be real.
851 for j
in xrange(i
+1):
852 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
854 Sij
= _embed_quaternion_matrix(Eij
)
857 # Beware, orthogonal but not normalized! The second,
858 # third, and fourth ones have a minus because they're
860 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
862 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
864 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
866 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
872 def _multiplication_table_from_matrix_basis(basis
):
874 At least three of the five simple Euclidean Jordan algebras have the
875 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
876 multiplication on the right is matrix multiplication. Given a basis
877 for the underlying matrix space, this function returns a
878 multiplication table (obtained by looping through the basis
879 elements) for an algebra of those matrices.
881 # In S^2, for example, we nominally have four coordinates even
882 # though the space is of dimension three only. The vector space V
883 # is supposed to hold the entire long vector, and the subspace W
884 # of V will be spanned by the vectors that arise from symmetric
885 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
886 field
= basis
[0].base_ring()
887 dimension
= basis
[0].nrows()
889 V
= VectorSpace(field
, dimension
**2)
890 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
892 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
895 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
896 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
901 def _embed_complex_matrix(M
):
903 Embed the n-by-n complex matrix ``M`` into the space of real
904 matrices of size 2n-by-2n via the map the sends each entry `z = a +
905 bi` to the block matrix ``[[a,b],[-b,a]]``.
909 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
913 sage: F = QuadraticField(-1,'i')
914 sage: x1 = F(4 - 2*i)
915 sage: x2 = F(1 + 2*i)
918 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
919 sage: _embed_complex_matrix(M)
928 Embedding is a homomorphism (isomorphism, in fact)::
930 sage: set_random_seed()
931 sage: n = ZZ.random_element(5)
932 sage: F = QuadraticField(-1, 'i')
933 sage: X = random_matrix(F, n)
934 sage: Y = random_matrix(F, n)
935 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
936 sage: expected = _embed_complex_matrix(X*Y)
937 sage: actual == expected
943 raise ValueError("the matrix 'M' must be square")
944 field
= M
.base_ring()
949 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
951 # We can drop the imaginaries here.
952 return matrix
.block(field
.base_ring(), n
, blocks
)
955 def _unembed_complex_matrix(M
):
957 The inverse of _embed_complex_matrix().
961 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
962 ....: _unembed_complex_matrix)
966 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
967 ....: [-2, 1, -4, 3],
968 ....: [ 9, 10, 11, 12],
969 ....: [-10, 9, -12, 11] ])
970 sage: _unembed_complex_matrix(A)
972 [ 10*i + 9 12*i + 11]
976 Unembedding is the inverse of embedding::
978 sage: set_random_seed()
979 sage: F = QuadraticField(-1, 'i')
980 sage: M = random_matrix(F, 3)
981 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
987 raise ValueError("the matrix 'M' must be square")
988 if not n
.mod(2).is_zero():
989 raise ValueError("the matrix 'M' must be a complex embedding")
991 F
= QuadraticField(-1, 'i')
994 # Go top-left to bottom-right (reading order), converting every
995 # 2-by-2 block we see to a single complex element.
997 for k
in xrange(n
/2):
998 for j
in xrange(n
/2):
999 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1000 if submat
[0,0] != submat
[1,1]:
1001 raise ValueError('bad on-diagonal submatrix')
1002 if submat
[0,1] != -submat
[1,0]:
1003 raise ValueError('bad off-diagonal submatrix')
1004 z
= submat
[0,0] + submat
[0,1]*i
1007 return matrix(F
, n
/2, elements
)
1010 def _embed_quaternion_matrix(M
):
1012 Embed the n-by-n quaternion matrix ``M`` into the space of real
1013 matrices of size 4n-by-4n by first sending each quaternion entry
1014 `z = a + bi + cj + dk` to the block-complex matrix
1015 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1020 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
1024 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1025 sage: i,j,k = Q.gens()
1026 sage: x = 1 + 2*i + 3*j + 4*k
1027 sage: M = matrix(Q, 1, [[x]])
1028 sage: _embed_quaternion_matrix(M)
1034 Embedding is a homomorphism (isomorphism, in fact)::
1036 sage: set_random_seed()
1037 sage: n = ZZ.random_element(5)
1038 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1039 sage: X = random_matrix(Q, n)
1040 sage: Y = random_matrix(Q, n)
1041 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1042 sage: expected = _embed_quaternion_matrix(X*Y)
1043 sage: actual == expected
1047 quaternions
= M
.base_ring()
1050 raise ValueError("the matrix 'M' must be square")
1052 F
= QuadraticField(-1, 'i')
1057 t
= z
.coefficient_tuple()
1062 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1063 [-c
+ d
*i
, a
- b
*i
]])
1064 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1066 # We should have real entries by now, so use the realest field
1067 # we've got for the return value.
1068 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1071 def _unembed_quaternion_matrix(M
):
1073 The inverse of _embed_quaternion_matrix().
1077 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1078 ....: _unembed_quaternion_matrix)
1082 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1083 ....: [-2, 1, -4, 3],
1084 ....: [-3, 4, 1, -2],
1085 ....: [-4, -3, 2, 1]])
1086 sage: _unembed_quaternion_matrix(M)
1087 [1 + 2*i + 3*j + 4*k]
1091 Unembedding is the inverse of embedding::
1093 sage: set_random_seed()
1094 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1095 sage: M = random_matrix(Q, 3)
1096 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1102 raise ValueError("the matrix 'M' must be square")
1103 if not n
.mod(4).is_zero():
1104 raise ValueError("the matrix 'M' must be a complex embedding")
1106 Q
= QuaternionAlgebra(QQ
,-1,-1)
1109 # Go top-left to bottom-right (reading order), converting every
1110 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1113 for l
in xrange(n
/4):
1114 for m
in xrange(n
/4):
1115 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1116 if submat
[0,0] != submat
[1,1].conjugate():
1117 raise ValueError('bad on-diagonal submatrix')
1118 if submat
[0,1] != -submat
[1,0].conjugate():
1119 raise ValueError('bad off-diagonal submatrix')
1120 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1121 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1124 return matrix(Q
, n
/4, elements
)
1127 # The usual inner product on R^n.
1129 return x
.to_vector().inner_product(y
.to_vector())
1131 # The inner product used for the real symmetric simple EJA.
1132 # We keep it as a separate function because e.g. the complex
1133 # algebra uses the same inner product, except divided by 2.
1134 def _matrix_ip(X
,Y
):
1135 X_mat
= X
.natural_representation()
1136 Y_mat
= Y
.natural_representation()
1137 return (X_mat
*Y_mat
).trace()
1140 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1142 The rank-n simple EJA consisting of real symmetric n-by-n
1143 matrices, the usual symmetric Jordan product, and the trace inner
1144 product. It has dimension `(n^2 + n)/2` over the reals.
1148 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1152 sage: J = RealSymmetricEJA(2)
1153 sage: e0, e1, e2 = J.gens()
1163 The dimension of this algebra is `(n^2 + n) / 2`::
1165 sage: set_random_seed()
1166 sage: n = ZZ.random_element(1,5)
1167 sage: J = RealSymmetricEJA(n)
1168 sage: J.dimension() == (n^2 + n)/2
1171 The Jordan multiplication is what we think it is::
1173 sage: set_random_seed()
1174 sage: n = ZZ.random_element(1,5)
1175 sage: J = RealSymmetricEJA(n)
1176 sage: x = J.random_element()
1177 sage: y = J.random_element()
1178 sage: actual = (x*y).natural_representation()
1179 sage: X = x.natural_representation()
1180 sage: Y = y.natural_representation()
1181 sage: expected = (X*Y + Y*X)/2
1182 sage: actual == expected
1184 sage: J(expected) == x*y
1187 We can change the generator prefix::
1189 sage: RealSymmetricEJA(3, prefix='q').gens()
1190 (q0, q1, q2, q3, q4, q5)
1193 def __init__(self
, n
, field
=QQ
, **kwargs
):
1194 S
= _real_symmetric_basis(n
, field
=field
)
1195 Qs
= _multiplication_table_from_matrix_basis(S
)
1197 fdeja
= super(RealSymmetricEJA
, self
)
1198 return fdeja
.__init
__(field
,
1204 def inner_product(self
, x
, y
):
1205 return _matrix_ip(x
,y
)
1208 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1210 The rank-n simple EJA consisting of complex Hermitian n-by-n
1211 matrices over the real numbers, the usual symmetric Jordan product,
1212 and the real-part-of-trace inner product. It has dimension `n^2` over
1217 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1221 The dimension of this algebra is `n^2`::
1223 sage: set_random_seed()
1224 sage: n = ZZ.random_element(1,5)
1225 sage: J = ComplexHermitianEJA(n)
1226 sage: J.dimension() == 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 = ComplexHermitianEJA(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: ComplexHermitianEJA(2, prefix='z').gens()
1251 def __init__(self
, n
, field
=QQ
, **kwargs
):
1252 S
= _complex_hermitian_basis(n
)
1253 Qs
= _multiplication_table_from_matrix_basis(S
)
1255 fdeja
= super(ComplexHermitianEJA
, self
)
1256 return fdeja
.__init
__(field
,
1263 def inner_product(self
, x
, y
):
1264 # Since a+bi on the diagonal is represented as
1269 # we'll double-count the "a" entries if we take the trace of
1271 return _matrix_ip(x
,y
)/2
1274 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1276 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1277 matrices, the usual symmetric Jordan product, and the
1278 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1283 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1287 The dimension of this algebra is `n^2`::
1289 sage: set_random_seed()
1290 sage: n = ZZ.random_element(1,5)
1291 sage: J = QuaternionHermitianEJA(n)
1292 sage: J.dimension() == 2*(n^2) - n
1295 The Jordan multiplication is what we think it is::
1297 sage: set_random_seed()
1298 sage: n = ZZ.random_element(1,5)
1299 sage: J = QuaternionHermitianEJA(n)
1300 sage: x = J.random_element()
1301 sage: y = J.random_element()
1302 sage: actual = (x*y).natural_representation()
1303 sage: X = x.natural_representation()
1304 sage: Y = y.natural_representation()
1305 sage: expected = (X*Y + Y*X)/2
1306 sage: actual == expected
1308 sage: J(expected) == x*y
1311 We can change the generator prefix::
1313 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1314 (a0, a1, a2, a3, a4, a5)
1317 def __init__(self
, n
, field
=QQ
, **kwargs
):
1318 S
= _quaternion_hermitian_basis(n
)
1319 Qs
= _multiplication_table_from_matrix_basis(S
)
1321 fdeja
= super(QuaternionHermitianEJA
, self
)
1322 return fdeja
.__init
__(field
,
1328 def inner_product(self
, x
, y
):
1329 # Since a+bi+cj+dk on the diagonal is represented as
1331 # a + bi +cj + dk = [ a b c d]
1336 # we'll quadruple-count the "a" entries if we take the trace of
1338 return _matrix_ip(x
,y
)/4
1341 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1343 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1344 with the usual inner product and jordan product ``x*y =
1345 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1350 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1354 This multiplication table can be verified by hand::
1356 sage: J = JordanSpinEJA(4)
1357 sage: e0,e1,e2,e3 = J.gens()
1373 We can change the generator prefix::
1375 sage: JordanSpinEJA(2, prefix='B').gens()
1379 def __init__(self
, n
, field
=QQ
, **kwargs
):
1380 V
= VectorSpace(field
, n
)
1381 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1391 z0
= x
.inner_product(y
)
1392 zbar
= y0
*xbar
+ x0
*ybar
1393 z
= V([z0
] + zbar
.list())
1394 mult_table
[i
][j
] = z
1396 # The rank of the spin algebra is two, unless we're in a
1397 # one-dimensional ambient space (because the rank is bounded by
1398 # the ambient dimension).
1399 fdeja
= super(JordanSpinEJA
, self
)
1400 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1402 def inner_product(self
, x
, y
):
1403 return _usual_ip(x
,y
)