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 FreeModule
, 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 # Don't use the parent vector space directly here in case this
211 # happens to be a subalgebra. In that case, we would be e.g.
212 # two-dimensional but span_of_basis() would expect three
214 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
215 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
216 V1
= V
.span_of_basis( basis
)
217 b
= (V1
.basis() + V1
.complement().basis())
218 return V
.span_of_basis(b
)
222 def _charpoly_coeff(self
, i
):
224 Return the coefficient polynomial "a_{i}" of this algebra's
225 general characteristic polynomial.
227 Having this be a separate cached method lets us compute and
228 store the trace/determinant (a_{r-1} and a_{0} respectively)
229 separate from the entire characteristic polynomial.
231 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
232 R
= A_of_x
.base_ring()
234 # Guaranteed by theory
237 # Danger: the in-place modification is done for performance
238 # reasons (reconstructing a matrix with huge polynomial
239 # entries is slow), but I don't know how cached_method works,
240 # so it's highly possible that we're modifying some global
241 # list variable by reference, here. In other words, you
242 # probably shouldn't call this method twice on the same
243 # algebra, at the same time, in two threads
244 Ai_orig
= A_of_x
.column(i
)
245 A_of_x
.set_column(i
,xr
)
246 numerator
= A_of_x
.det()
247 A_of_x
.set_column(i
,Ai_orig
)
249 # We're relying on the theory here to ensure that each a_i is
250 # indeed back in R, and the added negative signs are to make
251 # the whole charpoly expression sum to zero.
252 return R(-numerator
/detA
)
256 def _charpoly_matrix_system(self
):
258 Compute the matrix whose entries A_ij are polynomials in
259 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
260 corresponding to `x^r` and the determinent of the matrix A =
261 [A_ij]. In other words, all of the fixed (cachable) data needed
262 to compute the coefficients of the characteristic polynomial.
267 # Turn my vector space into a module so that "vectors" can
268 # have multivatiate polynomial entries.
269 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
270 R
= PolynomialRing(self
.base_ring(), names
)
272 # Using change_ring() on the parent's vector space doesn't work
273 # here because, in a subalgebra, that vector space has a basis
274 # and change_ring() tries to bring the basis along with it. And
275 # that doesn't work unless the new ring is a PID, which it usually
279 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
283 # And figure out the "left multiplication by x" matrix in
286 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
287 for i
in range(n
) ] # don't recompute these!
289 ek
= self
.monomial(k
).to_vector()
291 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
292 for i
in range(n
) ) )
293 Lx
= matrix
.column(R
, lmbx_cols
)
295 # Now we can compute powers of x "symbolically"
296 x_powers
= [self
.one().to_vector(), x
]
297 for d
in range(2, r
+1):
298 x_powers
.append( Lx
*(x_powers
[-1]) )
300 idmat
= matrix
.identity(R
, n
)
302 W
= self
._charpoly
_basis
_space
()
303 W
= W
.change_ring(R
.fraction_field())
305 # Starting with the standard coordinates x = (X1,X2,...,Xn)
306 # and then converting the entries to W-coordinates allows us
307 # to pass in the standard coordinates to the charpoly and get
308 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
311 # W.coordinates(x^2) eval'd at (standard z-coords)
315 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
317 # We want the middle equivalent thing in our matrix, but use
318 # the first equivalent thing instead so that we can pass in
319 # standard coordinates.
320 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
321 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
322 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
323 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
327 def characteristic_polynomial(self
):
329 Return a characteristic polynomial that works for all elements
332 The resulting polynomial has `n+1` variables, where `n` is the
333 dimension of this algebra. The first `n` variables correspond to
334 the coordinates of an algebra element: when evaluated at the
335 coordinates of an algebra element with respect to a certain
336 basis, the result is a univariate polynomial (in the one
337 remaining variable ``t``), namely the characteristic polynomial
342 sage: from mjo.eja.eja_algebra import JordanSpinEJA
346 The characteristic polynomial in the spin algebra is given in
347 Alizadeh, Example 11.11::
349 sage: J = JordanSpinEJA(3)
350 sage: p = J.characteristic_polynomial(); p
351 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
352 sage: xvec = J.one().to_vector()
360 # The list of coefficient polynomials a_1, a_2, ..., a_n.
361 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
363 # We go to a bit of trouble here to reorder the
364 # indeterminates, so that it's easier to evaluate the
365 # characteristic polynomial at x's coordinates and get back
366 # something in terms of t, which is what we want.
368 S
= PolynomialRing(self
.base_ring(),'t')
370 S
= PolynomialRing(S
, R
.variable_names())
373 # Note: all entries past the rth should be zero. The
374 # coefficient of the highest power (x^r) is 1, but it doesn't
375 # appear in the solution vector which contains coefficients
376 # for the other powers (to make them sum to x^r).
378 a
[r
] = 1 # corresponds to x^r
380 # When the rank is equal to the dimension, trying to
381 # assign a[r] goes out-of-bounds.
382 a
.append(1) # corresponds to x^r
384 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
387 def inner_product(self
, x
, y
):
389 The inner product associated with this Euclidean Jordan algebra.
391 Defaults to the trace inner product, but can be overridden by
392 subclasses if they are sure that the necessary properties are
397 sage: from mjo.eja.eja_algebra import random_eja
401 The inner product must satisfy its axiom for this algebra to truly
402 be a Euclidean Jordan Algebra::
404 sage: set_random_seed()
405 sage: J = random_eja()
406 sage: x = J.random_element()
407 sage: y = J.random_element()
408 sage: z = J.random_element()
409 sage: (x*y).inner_product(z) == y.inner_product(x*z)
413 if (not x
in self
) or (not y
in self
):
414 raise TypeError("arguments must live in this algebra")
415 return x
.trace_inner_product(y
)
418 def is_trivial(self
):
420 Return whether or not this algebra is trivial.
422 A trivial algebra contains only the zero element.
426 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
430 sage: J = ComplexHermitianEJA(3)
433 sage: A = J.zero().subalgebra_generated_by()
438 return self
.dimension() == 0
441 def multiplication_table(self
):
443 Return a visual representation of this algebra's multiplication
444 table (on basis elements).
448 sage: from mjo.eja.eja_algebra import JordanSpinEJA
452 sage: J = JordanSpinEJA(4)
453 sage: J.multiplication_table()
454 +----++----+----+----+----+
455 | * || e0 | e1 | e2 | e3 |
456 +====++====+====+====+====+
457 | e0 || e0 | e1 | e2 | e3 |
458 +----++----+----+----+----+
459 | e1 || e1 | e0 | 0 | 0 |
460 +----++----+----+----+----+
461 | e2 || e2 | 0 | e0 | 0 |
462 +----++----+----+----+----+
463 | e3 || e3 | 0 | 0 | e0 |
464 +----++----+----+----+----+
467 M
= list(self
._multiplication
_table
) # copy
468 for i
in range(len(M
)):
469 # M had better be "square"
470 M
[i
] = [self
.monomial(i
)] + M
[i
]
471 M
= [["*"] + list(self
.gens())] + M
472 return table(M
, header_row
=True, header_column
=True, frame
=True)
475 def natural_basis(self
):
477 Return a more-natural representation of this algebra's basis.
479 Every finite-dimensional Euclidean Jordan Algebra is a direct
480 sum of five simple algebras, four of which comprise Hermitian
481 matrices. This method returns the original "natural" basis
482 for our underlying vector space. (Typically, the natural basis
483 is used to construct the multiplication table in the first place.)
485 Note that this will always return a matrix. The standard basis
486 in `R^n` will be returned as `n`-by-`1` column matrices.
490 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
491 ....: RealSymmetricEJA)
495 sage: J = RealSymmetricEJA(2)
497 Finite family {0: e0, 1: e1, 2: e2}
498 sage: J.natural_basis()
506 sage: J = JordanSpinEJA(2)
508 Finite family {0: e0, 1: e1}
509 sage: J.natural_basis()
516 if self
._natural
_basis
is None:
517 return tuple( b
.to_vector().column() for b
in self
.basis() )
519 return self
._natural
_basis
525 Return the unit element of this algebra.
529 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
534 sage: J = RealCartesianProductEJA(5)
536 e0 + e1 + e2 + e3 + e4
540 The identity element acts like the identity::
542 sage: set_random_seed()
543 sage: J = random_eja()
544 sage: x = J.random_element()
545 sage: J.one()*x == x and x*J.one() == x
548 The matrix of the unit element's operator is the identity::
550 sage: set_random_seed()
551 sage: J = random_eja()
552 sage: actual = J.one().operator().matrix()
553 sage: expected = matrix.identity(J.base_ring(), J.dimension())
554 sage: actual == expected
558 # We can brute-force compute the matrices of the operators
559 # that correspond to the basis elements of this algebra.
560 # If some linear combination of those basis elements is the
561 # algebra identity, then the same linear combination of
562 # their matrices has to be the identity matrix.
564 # Of course, matrices aren't vectors in sage, so we have to
565 # appeal to the "long vectors" isometry.
566 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
568 # Now we use basis linear algebra to find the coefficients,
569 # of the matrices-as-vectors-linear-combination, which should
570 # work for the original algebra basis too.
571 A
= matrix
.column(self
.base_ring(), oper_vecs
)
573 # We used the isometry on the left-hand side already, but we
574 # still need to do it for the right-hand side. Recall that we
575 # wanted something that summed to the identity matrix.
576 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
578 # Now if there's an identity element in the algebra, this should work.
579 coeffs
= A
.solve_right(b
)
580 return self
.linear_combination(zip(self
.gens(), coeffs
))
583 def random_element(self
):
584 # Temporary workaround for https://trac.sagemath.org/ticket/28327
585 if self
.is_trivial():
588 s
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
589 return s
.random_element()
594 Return the rank of this EJA.
598 The author knows of no algorithm to compute the rank of an EJA
599 where only the multiplication table is known. In lieu of one, we
600 require the rank to be specified when the algebra is created,
601 and simply pass along that number here.
605 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
606 ....: RealSymmetricEJA,
607 ....: ComplexHermitianEJA,
608 ....: QuaternionHermitianEJA,
613 The rank of the Jordan spin algebra is always two::
615 sage: JordanSpinEJA(2).rank()
617 sage: JordanSpinEJA(3).rank()
619 sage: JordanSpinEJA(4).rank()
622 The rank of the `n`-by-`n` Hermitian real, complex, or
623 quaternion matrices is `n`::
625 sage: RealSymmetricEJA(2).rank()
627 sage: ComplexHermitianEJA(2).rank()
629 sage: QuaternionHermitianEJA(2).rank()
631 sage: RealSymmetricEJA(5).rank()
633 sage: ComplexHermitianEJA(5).rank()
635 sage: QuaternionHermitianEJA(5).rank()
640 Ensure that every EJA that we know how to construct has a
641 positive integer rank::
643 sage: set_random_seed()
644 sage: r = random_eja().rank()
645 sage: r in ZZ and r > 0
652 def vector_space(self
):
654 Return the vector space that underlies this algebra.
658 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
662 sage: J = RealSymmetricEJA(2)
663 sage: J.vector_space()
664 Vector space of dimension 3 over Rational Field
667 return self
.zero().to_vector().parent().ambient_vector_space()
670 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
673 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
675 Return the Euclidean Jordan Algebra corresponding to the set
676 `R^n` under the Hadamard product.
678 Note: this is nothing more than the Cartesian product of ``n``
679 copies of the spin algebra. Once Cartesian product algebras
680 are implemented, this can go.
684 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
688 This multiplication table can be verified by hand::
690 sage: J = RealCartesianProductEJA(3)
691 sage: e0,e1,e2 = J.gens()
707 We can change the generator prefix::
709 sage: RealCartesianProductEJA(3, prefix='r').gens()
713 def __init__(self
, n
, field
=QQ
, **kwargs
):
714 V
= VectorSpace(field
, n
)
715 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
718 fdeja
= super(RealCartesianProductEJA
, self
)
719 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
721 def inner_product(self
, x
, y
):
722 return _usual_ip(x
,y
)
727 Return a "random" finite-dimensional Euclidean Jordan Algebra.
731 For now, we choose a random natural number ``n`` (greater than zero)
732 and then give you back one of the following:
734 * The cartesian product of the rational numbers ``n`` times; this is
735 ``QQ^n`` with the Hadamard product.
737 * The Jordan spin algebra on ``QQ^n``.
739 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
742 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
743 in the space of ``2n``-by-``2n`` real symmetric matrices.
745 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
746 in the space of ``4n``-by-``4n`` real symmetric matrices.
748 Later this might be extended to return Cartesian products of the
753 sage: from mjo.eja.eja_algebra import random_eja
758 Euclidean Jordan algebra of dimension...
762 # The max_n component lets us choose different upper bounds on the
763 # value "n" that gets passed to the constructor. This is needed
764 # because e.g. R^{10} is reasonable to test, while the Hermitian
765 # 10-by-10 quaternion matrices are not.
766 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
768 (RealSymmetricEJA
, 5),
769 (ComplexHermitianEJA
, 4),
770 (QuaternionHermitianEJA
, 3)])
771 n
= ZZ
.random_element(1, max_n
)
772 return constructor(n
, field
=QQ
)
776 def _real_symmetric_basis(n
, field
=QQ
):
778 Return a basis for the space of real symmetric n-by-n matrices.
780 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
784 for j
in xrange(i
+1):
785 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
789 # Beware, orthogonal but not normalized!
790 Sij
= Eij
+ Eij
.transpose()
795 def _complex_hermitian_basis(n
, field
=QQ
):
797 Returns a basis for the space of complex Hermitian n-by-n matrices.
801 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
805 sage: set_random_seed()
806 sage: n = ZZ.random_element(1,5)
807 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
811 F
= QuadraticField(-1, 'I')
814 # This is like the symmetric case, but we need to be careful:
816 # * We want conjugate-symmetry, not just symmetry.
817 # * The diagonal will (as a result) be real.
821 for j
in xrange(i
+1):
822 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
824 Sij
= _embed_complex_matrix(Eij
)
827 # Beware, orthogonal but not normalized! The second one
828 # has a minus because it's conjugated.
829 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
831 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
836 def _quaternion_hermitian_basis(n
, field
=QQ
):
838 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
842 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
846 sage: set_random_seed()
847 sage: n = ZZ.random_element(1,5)
848 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
852 Q
= QuaternionAlgebra(QQ
,-1,-1)
855 # This is like the symmetric case, but we need to be careful:
857 # * We want conjugate-symmetry, not just symmetry.
858 # * The diagonal will (as a result) be real.
862 for j
in xrange(i
+1):
863 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
865 Sij
= _embed_quaternion_matrix(Eij
)
868 # Beware, orthogonal but not normalized! The second,
869 # third, and fourth ones have a minus because they're
871 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
873 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
875 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
877 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
883 def _multiplication_table_from_matrix_basis(basis
):
885 At least three of the five simple Euclidean Jordan algebras have the
886 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
887 multiplication on the right is matrix multiplication. Given a basis
888 for the underlying matrix space, this function returns a
889 multiplication table (obtained by looping through the basis
890 elements) for an algebra of those matrices.
892 # In S^2, for example, we nominally have four coordinates even
893 # though the space is of dimension three only. The vector space V
894 # is supposed to hold the entire long vector, and the subspace W
895 # of V will be spanned by the vectors that arise from symmetric
896 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
897 field
= basis
[0].base_ring()
898 dimension
= basis
[0].nrows()
900 V
= VectorSpace(field
, dimension
**2)
901 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
903 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
906 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
907 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
912 def _embed_complex_matrix(M
):
914 Embed the n-by-n complex matrix ``M`` into the space of real
915 matrices of size 2n-by-2n via the map the sends each entry `z = a +
916 bi` to the block matrix ``[[a,b],[-b,a]]``.
920 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
924 sage: F = QuadraticField(-1,'i')
925 sage: x1 = F(4 - 2*i)
926 sage: x2 = F(1 + 2*i)
929 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
930 sage: _embed_complex_matrix(M)
939 Embedding is a homomorphism (isomorphism, in fact)::
941 sage: set_random_seed()
942 sage: n = ZZ.random_element(5)
943 sage: F = QuadraticField(-1, 'i')
944 sage: X = random_matrix(F, n)
945 sage: Y = random_matrix(F, n)
946 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
947 sage: expected = _embed_complex_matrix(X*Y)
948 sage: actual == expected
954 raise ValueError("the matrix 'M' must be square")
955 field
= M
.base_ring()
960 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
962 # We can drop the imaginaries here.
963 return matrix
.block(field
.base_ring(), n
, blocks
)
966 def _unembed_complex_matrix(M
):
968 The inverse of _embed_complex_matrix().
972 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
973 ....: _unembed_complex_matrix)
977 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
978 ....: [-2, 1, -4, 3],
979 ....: [ 9, 10, 11, 12],
980 ....: [-10, 9, -12, 11] ])
981 sage: _unembed_complex_matrix(A)
983 [ 10*i + 9 12*i + 11]
987 Unembedding is the inverse of embedding::
989 sage: set_random_seed()
990 sage: F = QuadraticField(-1, 'i')
991 sage: M = random_matrix(F, 3)
992 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
998 raise ValueError("the matrix 'M' must be square")
999 if not n
.mod(2).is_zero():
1000 raise ValueError("the matrix 'M' must be a complex embedding")
1002 F
= QuadraticField(-1, 'i')
1005 # Go top-left to bottom-right (reading order), converting every
1006 # 2-by-2 block we see to a single complex element.
1008 for k
in xrange(n
/2):
1009 for j
in xrange(n
/2):
1010 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1011 if submat
[0,0] != submat
[1,1]:
1012 raise ValueError('bad on-diagonal submatrix')
1013 if submat
[0,1] != -submat
[1,0]:
1014 raise ValueError('bad off-diagonal submatrix')
1015 z
= submat
[0,0] + submat
[0,1]*i
1018 return matrix(F
, n
/2, elements
)
1021 def _embed_quaternion_matrix(M
):
1023 Embed the n-by-n quaternion matrix ``M`` into the space of real
1024 matrices of size 4n-by-4n by first sending each quaternion entry
1025 `z = a + bi + cj + dk` to the block-complex matrix
1026 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1031 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
1035 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1036 sage: i,j,k = Q.gens()
1037 sage: x = 1 + 2*i + 3*j + 4*k
1038 sage: M = matrix(Q, 1, [[x]])
1039 sage: _embed_quaternion_matrix(M)
1045 Embedding is a homomorphism (isomorphism, in fact)::
1047 sage: set_random_seed()
1048 sage: n = ZZ.random_element(5)
1049 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1050 sage: X = random_matrix(Q, n)
1051 sage: Y = random_matrix(Q, n)
1052 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1053 sage: expected = _embed_quaternion_matrix(X*Y)
1054 sage: actual == expected
1058 quaternions
= M
.base_ring()
1061 raise ValueError("the matrix 'M' must be square")
1063 F
= QuadraticField(-1, 'i')
1068 t
= z
.coefficient_tuple()
1073 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1074 [-c
+ d
*i
, a
- b
*i
]])
1075 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1077 # We should have real entries by now, so use the realest field
1078 # we've got for the return value.
1079 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1082 def _unembed_quaternion_matrix(M
):
1084 The inverse of _embed_quaternion_matrix().
1088 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1089 ....: _unembed_quaternion_matrix)
1093 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1094 ....: [-2, 1, -4, 3],
1095 ....: [-3, 4, 1, -2],
1096 ....: [-4, -3, 2, 1]])
1097 sage: _unembed_quaternion_matrix(M)
1098 [1 + 2*i + 3*j + 4*k]
1102 Unembedding is the inverse of embedding::
1104 sage: set_random_seed()
1105 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1106 sage: M = random_matrix(Q, 3)
1107 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1113 raise ValueError("the matrix 'M' must be square")
1114 if not n
.mod(4).is_zero():
1115 raise ValueError("the matrix 'M' must be a complex embedding")
1117 Q
= QuaternionAlgebra(QQ
,-1,-1)
1120 # Go top-left to bottom-right (reading order), converting every
1121 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1124 for l
in xrange(n
/4):
1125 for m
in xrange(n
/4):
1126 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1127 if submat
[0,0] != submat
[1,1].conjugate():
1128 raise ValueError('bad on-diagonal submatrix')
1129 if submat
[0,1] != -submat
[1,0].conjugate():
1130 raise ValueError('bad off-diagonal submatrix')
1131 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1132 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1135 return matrix(Q
, n
/4, elements
)
1138 # The usual inner product on R^n.
1140 return x
.to_vector().inner_product(y
.to_vector())
1142 # The inner product used for the real symmetric simple EJA.
1143 # We keep it as a separate function because e.g. the complex
1144 # algebra uses the same inner product, except divided by 2.
1145 def _matrix_ip(X
,Y
):
1146 X_mat
= X
.natural_representation()
1147 Y_mat
= Y
.natural_representation()
1148 return (X_mat
*Y_mat
).trace()
1151 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1153 The rank-n simple EJA consisting of real symmetric n-by-n
1154 matrices, the usual symmetric Jordan product, and the trace inner
1155 product. It has dimension `(n^2 + n)/2` over the reals.
1159 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1163 sage: J = RealSymmetricEJA(2)
1164 sage: e0, e1, e2 = J.gens()
1174 The dimension of this algebra is `(n^2 + n) / 2`::
1176 sage: set_random_seed()
1177 sage: n = ZZ.random_element(1,5)
1178 sage: J = RealSymmetricEJA(n)
1179 sage: J.dimension() == (n^2 + n)/2
1182 The Jordan multiplication is what we think it is::
1184 sage: set_random_seed()
1185 sage: n = ZZ.random_element(1,5)
1186 sage: J = RealSymmetricEJA(n)
1187 sage: x = J.random_element()
1188 sage: y = J.random_element()
1189 sage: actual = (x*y).natural_representation()
1190 sage: X = x.natural_representation()
1191 sage: Y = y.natural_representation()
1192 sage: expected = (X*Y + Y*X)/2
1193 sage: actual == expected
1195 sage: J(expected) == x*y
1198 We can change the generator prefix::
1200 sage: RealSymmetricEJA(3, prefix='q').gens()
1201 (q0, q1, q2, q3, q4, q5)
1204 def __init__(self
, n
, field
=QQ
, **kwargs
):
1205 S
= _real_symmetric_basis(n
, field
=field
)
1206 Qs
= _multiplication_table_from_matrix_basis(S
)
1208 fdeja
= super(RealSymmetricEJA
, self
)
1209 return fdeja
.__init
__(field
,
1215 def inner_product(self
, x
, y
):
1216 return _matrix_ip(x
,y
)
1219 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1221 The rank-n simple EJA consisting of complex Hermitian n-by-n
1222 matrices over the real numbers, the usual symmetric Jordan product,
1223 and the real-part-of-trace inner product. It has dimension `n^2` over
1228 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1232 The dimension of this algebra is `n^2`::
1234 sage: set_random_seed()
1235 sage: n = ZZ.random_element(1,5)
1236 sage: J = ComplexHermitianEJA(n)
1237 sage: J.dimension() == n^2
1240 The Jordan multiplication is what we think it is::
1242 sage: set_random_seed()
1243 sage: n = ZZ.random_element(1,5)
1244 sage: J = ComplexHermitianEJA(n)
1245 sage: x = J.random_element()
1246 sage: y = J.random_element()
1247 sage: actual = (x*y).natural_representation()
1248 sage: X = x.natural_representation()
1249 sage: Y = y.natural_representation()
1250 sage: expected = (X*Y + Y*X)/2
1251 sage: actual == expected
1253 sage: J(expected) == x*y
1256 We can change the generator prefix::
1258 sage: ComplexHermitianEJA(2, prefix='z').gens()
1262 def __init__(self
, n
, field
=QQ
, **kwargs
):
1263 S
= _complex_hermitian_basis(n
)
1264 Qs
= _multiplication_table_from_matrix_basis(S
)
1266 fdeja
= super(ComplexHermitianEJA
, self
)
1267 return fdeja
.__init
__(field
,
1274 def inner_product(self
, x
, y
):
1275 # Since a+bi on the diagonal is represented as
1280 # we'll double-count the "a" entries if we take the trace of
1282 return _matrix_ip(x
,y
)/2
1285 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1287 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1288 matrices, the usual symmetric Jordan product, and the
1289 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1294 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1298 The dimension of this algebra is `n^2`::
1300 sage: set_random_seed()
1301 sage: n = ZZ.random_element(1,5)
1302 sage: J = QuaternionHermitianEJA(n)
1303 sage: J.dimension() == 2*(n^2) - n
1306 The Jordan multiplication is what we think it is::
1308 sage: set_random_seed()
1309 sage: n = ZZ.random_element(1,5)
1310 sage: J = QuaternionHermitianEJA(n)
1311 sage: x = J.random_element()
1312 sage: y = J.random_element()
1313 sage: actual = (x*y).natural_representation()
1314 sage: X = x.natural_representation()
1315 sage: Y = y.natural_representation()
1316 sage: expected = (X*Y + Y*X)/2
1317 sage: actual == expected
1319 sage: J(expected) == x*y
1322 We can change the generator prefix::
1324 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1325 (a0, a1, a2, a3, a4, a5)
1328 def __init__(self
, n
, field
=QQ
, **kwargs
):
1329 S
= _quaternion_hermitian_basis(n
)
1330 Qs
= _multiplication_table_from_matrix_basis(S
)
1332 fdeja
= super(QuaternionHermitianEJA
, self
)
1333 return fdeja
.__init
__(field
,
1339 def inner_product(self
, x
, y
):
1340 # Since a+bi+cj+dk on the diagonal is represented as
1342 # a + bi +cj + dk = [ a b c d]
1347 # we'll quadruple-count the "a" entries if we take the trace of
1349 return _matrix_ip(x
,y
)/4
1352 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1354 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1355 with the usual inner product and jordan product ``x*y =
1356 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1361 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1365 This multiplication table can be verified by hand::
1367 sage: J = JordanSpinEJA(4)
1368 sage: e0,e1,e2,e3 = J.gens()
1384 We can change the generator prefix::
1386 sage: JordanSpinEJA(2, prefix='B').gens()
1390 def __init__(self
, n
, field
=QQ
, **kwargs
):
1391 V
= VectorSpace(field
, n
)
1392 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1402 z0
= x
.inner_product(y
)
1403 zbar
= y0
*xbar
+ x0
*ybar
1404 z
= V([z0
] + zbar
.list())
1405 mult_table
[i
][j
] = z
1407 # The rank of the spin algebra is two, unless we're in a
1408 # one-dimensional ambient space (because the rank is bounded by
1409 # the ambient dimension).
1410 fdeja
= super(JordanSpinEJA
, self
)
1411 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1413 def inner_product(self
, x
, y
):
1414 return _usual_ip(x
,y
)