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 itertools
import repeat
10 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
11 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
12 from sage
.combinat
.free_module
import CombinatorialFreeModule
13 from sage
.matrix
.constructor
import matrix
14 from sage
.matrix
.matrix_space
import MatrixSpace
15 from sage
.misc
.cachefunc
import cached_method
16 from sage
.misc
.lazy_import
import lazy_import
17 from sage
.misc
.prandom
import choice
18 from sage
.misc
.table
import table
19 from sage
.modules
.free_module
import FreeModule
, VectorSpace
20 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
23 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
24 lazy_import('mjo.eja.eja_subalgebra',
25 'FiniteDimensionalEuclideanJordanSubalgebra')
26 from mjo
.eja
.eja_utils
import _mat2vec
28 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
30 def _coerce_map_from_base_ring(self
):
32 Disable the map from the base ring into the algebra.
34 Performing a nonsense conversion like this automatically
35 is counterpedagogical. The fallback is to try the usual
36 element constructor, which should also fail.
40 sage: from mjo.eja.eja_algebra import random_eja
44 sage: set_random_seed()
45 sage: J = random_eja()
47 Traceback (most recent call last):
49 ValueError: not a naturally-represented algebra element
64 sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
68 By definition, Jordan multiplication commutes::
70 sage: set_random_seed()
71 sage: J = random_eja()
72 sage: x,y = J.random_elements(2)
78 The ``field`` we're given must be real::
80 sage: JordanSpinEJA(2,QQbar)
81 Traceback (most recent call last):
83 ValueError: field is not real
87 if not field
.is_subring(RR
):
88 # Note: this does return true for the real algebraic
89 # field, and any quadratic field where we've specified
91 raise ValueError('field is not real')
93 self
._natural
_basis
= natural_basis
96 category
= MagmaticAlgebras(field
).FiniteDimensional()
97 category
= category
.WithBasis().Unital()
99 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
101 range(len(mult_table
)),
104 self
.print_options(bracket
='')
106 # The multiplication table we're given is necessarily in terms
107 # of vectors, because we don't have an algebra yet for
108 # anything to be an element of. However, it's faster in the
109 # long run to have the multiplication table be in terms of
110 # algebra elements. We do this after calling the superclass
111 # constructor so that from_vector() knows what to do.
112 self
._multiplication
_table
= [
113 list(map(lambda x
: self
.from_vector(x
), ls
))
118 def _element_constructor_(self
, elt
):
120 Construct an element of this algebra from its natural
123 This gets called only after the parent element _call_ method
124 fails to find a coercion for the argument.
128 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
130 ....: RealSymmetricEJA)
134 The identity in `S^n` is converted to the identity in the EJA::
136 sage: J = RealSymmetricEJA(3)
137 sage: I = matrix.identity(QQ,3)
138 sage: J(I) == J.one()
141 This skew-symmetric matrix can't be represented in the EJA::
143 sage: J = RealSymmetricEJA(3)
144 sage: A = matrix(QQ,3, lambda i,j: i-j)
146 Traceback (most recent call last):
148 ArithmeticError: vector is not in free module
152 Ensure that we can convert any element of the two non-matrix
153 simple algebras (whose natural representations are their usual
154 vector representations) back and forth faithfully::
156 sage: set_random_seed()
157 sage: J = HadamardEJA.random_instance()
158 sage: x = J.random_element()
159 sage: J(x.to_vector().column()) == x
161 sage: J = JordanSpinEJA.random_instance()
162 sage: x = J.random_element()
163 sage: J(x.to_vector().column()) == x
167 msg
= "not a naturally-represented algebra element"
169 # The superclass implementation of random_element()
170 # needs to be able to coerce "0" into the algebra.
172 elif elt
in self
.base_ring():
173 # Ensure that no base ring -> algebra coercion is performed
174 # by this method. There's some stupidity in sage that would
175 # otherwise propagate to this method; for example, sage thinks
176 # that the integer 3 belongs to the space of 2-by-2 matrices.
177 raise ValueError(msg
)
179 natural_basis
= self
.natural_basis()
180 basis_space
= natural_basis
[0].matrix_space()
181 if elt
not in basis_space
:
182 raise ValueError(msg
)
184 # Thanks for nothing! Matrix spaces aren't vector spaces in
185 # Sage, so we have to figure out its natural-basis coordinates
186 # ourselves. We use the basis space's ring instead of the
187 # element's ring because the basis space might be an algebraic
188 # closure whereas the base ring of the 3-by-3 identity matrix
189 # could be QQ instead of QQbar.
190 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
191 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
192 coords
= W
.coordinate_vector(_mat2vec(elt
))
193 return self
.from_vector(coords
)
198 Return a string representation of ``self``.
202 sage: from mjo.eja.eja_algebra import JordanSpinEJA
206 Ensure that it says what we think it says::
208 sage: JordanSpinEJA(2, field=AA)
209 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
210 sage: JordanSpinEJA(3, field=RDF)
211 Euclidean Jordan algebra of dimension 3 over Real Double Field
214 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
215 return fmt
.format(self
.dimension(), self
.base_ring())
217 def product_on_basis(self
, i
, j
):
218 return self
._multiplication
_table
[i
][j
]
220 def _a_regular_element(self
):
222 Guess a regular element. Needed to compute the basis for our
223 characteristic polynomial coefficients.
227 sage: from mjo.eja.eja_algebra import random_eja
231 Ensure that this hacky method succeeds for every algebra that we
232 know how to construct::
234 sage: set_random_seed()
235 sage: J = random_eja()
236 sage: J._a_regular_element().is_regular()
241 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
242 if not z
.is_regular():
243 raise ValueError("don't know a regular element")
248 def _charpoly_basis_space(self
):
250 Return the vector space spanned by the basis used in our
251 characteristic polynomial coefficients. This is used not only to
252 compute those coefficients, but also any time we need to
253 evaluate the coefficients (like when we compute the trace or
256 z
= self
._a
_regular
_element
()
257 # Don't use the parent vector space directly here in case this
258 # happens to be a subalgebra. In that case, we would be e.g.
259 # two-dimensional but span_of_basis() would expect three
261 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
262 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
263 V1
= V
.span_of_basis( basis
)
264 b
= (V1
.basis() + V1
.complement().basis())
265 return V
.span_of_basis(b
)
270 def _charpoly_coeff(self
, i
):
272 Return the coefficient polynomial "a_{i}" of this algebra's
273 general characteristic polynomial.
275 Having this be a separate cached method lets us compute and
276 store the trace/determinant (a_{r-1} and a_{0} respectively)
277 separate from the entire characteristic polynomial.
279 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
280 R
= A_of_x
.base_ring()
285 # Guaranteed by theory
288 # Danger: the in-place modification is done for performance
289 # reasons (reconstructing a matrix with huge polynomial
290 # entries is slow), but I don't know how cached_method works,
291 # so it's highly possible that we're modifying some global
292 # list variable by reference, here. In other words, you
293 # probably shouldn't call this method twice on the same
294 # algebra, at the same time, in two threads
295 Ai_orig
= A_of_x
.column(i
)
296 A_of_x
.set_column(i
,xr
)
297 numerator
= A_of_x
.det()
298 A_of_x
.set_column(i
,Ai_orig
)
300 # We're relying on the theory here to ensure that each a_i is
301 # indeed back in R, and the added negative signs are to make
302 # the whole charpoly expression sum to zero.
303 return R(-numerator
/detA
)
307 def _charpoly_matrix_system(self
):
309 Compute the matrix whose entries A_ij are polynomials in
310 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
311 corresponding to `x^r` and the determinent of the matrix A =
312 [A_ij]. In other words, all of the fixed (cachable) data needed
313 to compute the coefficients of the characteristic polynomial.
318 # Turn my vector space into a module so that "vectors" can
319 # have multivatiate polynomial entries.
320 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
321 R
= PolynomialRing(self
.base_ring(), names
)
323 # Using change_ring() on the parent's vector space doesn't work
324 # here because, in a subalgebra, that vector space has a basis
325 # and change_ring() tries to bring the basis along with it. And
326 # that doesn't work unless the new ring is a PID, which it usually
330 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
334 # And figure out the "left multiplication by x" matrix in
337 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
338 for i
in range(n
) ] # don't recompute these!
340 ek
= self
.monomial(k
).to_vector()
342 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
343 for i
in range(n
) ) )
344 Lx
= matrix
.column(R
, lmbx_cols
)
346 # Now we can compute powers of x "symbolically"
347 x_powers
= [self
.one().to_vector(), x
]
348 for d
in range(2, r
+1):
349 x_powers
.append( Lx
*(x_powers
[-1]) )
351 idmat
= matrix
.identity(R
, n
)
353 W
= self
._charpoly
_basis
_space
()
354 W
= W
.change_ring(R
.fraction_field())
356 # Starting with the standard coordinates x = (X1,X2,...,Xn)
357 # and then converting the entries to W-coordinates allows us
358 # to pass in the standard coordinates to the charpoly and get
359 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
362 # W.coordinates(x^2) eval'd at (standard z-coords)
366 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
368 # We want the middle equivalent thing in our matrix, but use
369 # the first equivalent thing instead so that we can pass in
370 # standard coordinates.
371 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
372 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
373 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
374 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
378 def characteristic_polynomial(self
):
380 Return a characteristic polynomial that works for all elements
383 The resulting polynomial has `n+1` variables, where `n` is the
384 dimension of this algebra. The first `n` variables correspond to
385 the coordinates of an algebra element: when evaluated at the
386 coordinates of an algebra element with respect to a certain
387 basis, the result is a univariate polynomial (in the one
388 remaining variable ``t``), namely the characteristic polynomial
393 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
397 The characteristic polynomial in the spin algebra is given in
398 Alizadeh, Example 11.11::
400 sage: J = JordanSpinEJA(3)
401 sage: p = J.characteristic_polynomial(); p
402 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
403 sage: xvec = J.one().to_vector()
407 By definition, the characteristic polynomial is a monic
408 degree-zero polynomial in a rank-zero algebra. Note that
409 Cayley-Hamilton is indeed satisfied since the polynomial
410 ``1`` evaluates to the identity element of the algebra on
413 sage: J = TrivialEJA()
414 sage: J.characteristic_polynomial()
421 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
422 a
= [ self
._charpoly
_coeff
(i
) for i
in range(r
+1) ]
424 # We go to a bit of trouble here to reorder the
425 # indeterminates, so that it's easier to evaluate the
426 # characteristic polynomial at x's coordinates and get back
427 # something in terms of t, which is what we want.
429 S
= PolynomialRing(self
.base_ring(),'t')
431 S
= PolynomialRing(S
, R
.variable_names())
434 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
437 def inner_product(self
, x
, y
):
439 The inner product associated with this Euclidean Jordan algebra.
441 Defaults to the trace inner product, but can be overridden by
442 subclasses if they are sure that the necessary properties are
447 sage: from mjo.eja.eja_algebra import random_eja
451 Our inner product is "associative," which means the following for
452 a symmetric bilinear form::
454 sage: set_random_seed()
455 sage: J = random_eja()
456 sage: x,y,z = J.random_elements(3)
457 sage: (x*y).inner_product(z) == y.inner_product(x*z)
461 X
= x
.natural_representation()
462 Y
= y
.natural_representation()
463 return self
.natural_inner_product(X
,Y
)
466 def is_trivial(self
):
468 Return whether or not this algebra is trivial.
470 A trivial algebra contains only the zero element.
474 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
479 sage: J = ComplexHermitianEJA(3)
485 sage: J = TrivialEJA()
490 return self
.dimension() == 0
493 def multiplication_table(self
):
495 Return a visual representation of this algebra's multiplication
496 table (on basis elements).
500 sage: from mjo.eja.eja_algebra import JordanSpinEJA
504 sage: J = JordanSpinEJA(4)
505 sage: J.multiplication_table()
506 +----++----+----+----+----+
507 | * || e0 | e1 | e2 | e3 |
508 +====++====+====+====+====+
509 | e0 || e0 | e1 | e2 | e3 |
510 +----++----+----+----+----+
511 | e1 || e1 | e0 | 0 | 0 |
512 +----++----+----+----+----+
513 | e2 || e2 | 0 | e0 | 0 |
514 +----++----+----+----+----+
515 | e3 || e3 | 0 | 0 | e0 |
516 +----++----+----+----+----+
519 M
= list(self
._multiplication
_table
) # copy
520 for i
in range(len(M
)):
521 # M had better be "square"
522 M
[i
] = [self
.monomial(i
)] + M
[i
]
523 M
= [["*"] + list(self
.gens())] + M
524 return table(M
, header_row
=True, header_column
=True, frame
=True)
527 def natural_basis(self
):
529 Return a more-natural representation of this algebra's basis.
531 Every finite-dimensional Euclidean Jordan Algebra is a direct
532 sum of five simple algebras, four of which comprise Hermitian
533 matrices. This method returns the original "natural" basis
534 for our underlying vector space. (Typically, the natural basis
535 is used to construct the multiplication table in the first place.)
537 Note that this will always return a matrix. The standard basis
538 in `R^n` will be returned as `n`-by-`1` column matrices.
542 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
543 ....: RealSymmetricEJA)
547 sage: J = RealSymmetricEJA(2)
549 Finite family {0: e0, 1: e1, 2: e2}
550 sage: J.natural_basis()
552 [1 0] [ 0 0.7071067811865475?] [0 0]
553 [0 0], [0.7071067811865475? 0], [0 1]
558 sage: J = JordanSpinEJA(2)
560 Finite family {0: e0, 1: e1}
561 sage: J.natural_basis()
568 if self
._natural
_basis
is None:
569 M
= self
.natural_basis_space()
570 return tuple( M(b
.to_vector()) for b
in self
.basis() )
572 return self
._natural
_basis
575 def natural_basis_space(self
):
577 Return the matrix space in which this algebra's natural basis
580 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
581 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
583 return self
._natural
_basis
[0].matrix_space()
587 def natural_inner_product(X
,Y
):
589 Compute the inner product of two naturally-represented elements.
591 For example in the real symmetric matrix EJA, this will compute
592 the trace inner-product of two n-by-n symmetric matrices. The
593 default should work for the real cartesian product EJA, the
594 Jordan spin EJA, and the real symmetric matrices. The others
595 will have to be overridden.
597 return (X
.conjugate_transpose()*Y
).trace()
603 Return the unit element of this algebra.
607 sage: from mjo.eja.eja_algebra import (HadamardEJA,
612 sage: J = HadamardEJA(5)
614 e0 + e1 + e2 + e3 + e4
618 The identity element acts like the identity::
620 sage: set_random_seed()
621 sage: J = random_eja()
622 sage: x = J.random_element()
623 sage: J.one()*x == x and x*J.one() == x
626 The matrix of the unit element's operator is the identity::
628 sage: set_random_seed()
629 sage: J = random_eja()
630 sage: actual = J.one().operator().matrix()
631 sage: expected = matrix.identity(J.base_ring(), J.dimension())
632 sage: actual == expected
636 # We can brute-force compute the matrices of the operators
637 # that correspond to the basis elements of this algebra.
638 # If some linear combination of those basis elements is the
639 # algebra identity, then the same linear combination of
640 # their matrices has to be the identity matrix.
642 # Of course, matrices aren't vectors in sage, so we have to
643 # appeal to the "long vectors" isometry.
644 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
646 # Now we use basis linear algebra to find the coefficients,
647 # of the matrices-as-vectors-linear-combination, which should
648 # work for the original algebra basis too.
649 A
= matrix
.column(self
.base_ring(), oper_vecs
)
651 # We used the isometry on the left-hand side already, but we
652 # still need to do it for the right-hand side. Recall that we
653 # wanted something that summed to the identity matrix.
654 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
656 # Now if there's an identity element in the algebra, this should work.
657 coeffs
= A
.solve_right(b
)
658 return self
.linear_combination(zip(self
.gens(), coeffs
))
661 def peirce_decomposition(self
, c
):
663 The Peirce decomposition of this algebra relative to the
666 In the future, this can be extended to a complete system of
667 orthogonal idempotents.
671 - ``c`` -- an idempotent of this algebra.
675 A triple (J0, J5, J1) containing two subalgebras and one subspace
678 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
679 corresponding to the eigenvalue zero.
681 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
682 corresponding to the eigenvalue one-half.
684 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
685 corresponding to the eigenvalue one.
687 These are the only possible eigenspaces for that operator, and this
688 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
689 orthogonal, and are subalgebras of this algebra with the appropriate
694 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
698 The canonical example comes from the symmetric matrices, which
699 decompose into diagonal and off-diagonal parts::
701 sage: J = RealSymmetricEJA(3)
702 sage: C = matrix(QQ, [ [1,0,0],
706 sage: J0,J5,J1 = J.peirce_decomposition(c)
708 Euclidean Jordan algebra of dimension 1...
710 Vector space of degree 6 and dimension 2...
712 Euclidean Jordan algebra of dimension 3...
716 Every algebra decomposes trivially with respect to its identity
719 sage: set_random_seed()
720 sage: J = random_eja()
721 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
722 sage: J0.dimension() == 0 and J5.dimension() == 0
724 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
727 The identity elements in the two subalgebras are the
728 projections onto their respective subspaces of the
729 superalgebra's identity element::
731 sage: set_random_seed()
732 sage: J = random_eja()
733 sage: x = J.random_element()
734 sage: if not J.is_trivial():
735 ....: while x.is_nilpotent():
736 ....: x = J.random_element()
737 sage: c = x.subalgebra_idempotent()
738 sage: J0,J5,J1 = J.peirce_decomposition(c)
739 sage: J1(c) == J1.one()
741 sage: J0(J.one() - c) == J0.one()
745 if not c
.is_idempotent():
746 raise ValueError("element is not idempotent: %s" % c
)
748 # Default these to what they should be if they turn out to be
749 # trivial, because eigenspaces_left() won't return eigenvalues
750 # corresponding to trivial spaces (e.g. it returns only the
751 # eigenspace corresponding to lambda=1 if you take the
752 # decomposition relative to the identity element).
753 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
754 J0
= trivial
# eigenvalue zero
755 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
756 J1
= trivial
# eigenvalue one
758 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
759 if eigval
== ~
(self
.base_ring()(2)):
762 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
763 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
769 raise ValueError("unexpected eigenvalue: %s" % eigval
)
774 def random_elements(self
, count
):
776 Return ``count`` random elements as a tuple.
780 sage: from mjo.eja.eja_algebra import JordanSpinEJA
784 sage: J = JordanSpinEJA(3)
785 sage: x,y,z = J.random_elements(3)
786 sage: all( [ x in J, y in J, z in J ])
788 sage: len( J.random_elements(10) ) == 10
792 return tuple( self
.random_element() for idx
in range(count
) )
797 Return the rank of this EJA.
801 We first compute the polynomial "column matrices" `p_{k}` that
802 evaluate to `x^k` on the coordinates of `x`. Then, we begin
803 adding them to a matrix one at a time, and trying to solve the
804 system that makes `p_{0}`,`p_{1}`,..., `p_{s-1}` add up to
805 `p_{s}`. This will succeed only when `s` is the rank of the
806 algebra, as proven in a recent draft paper of mine.
810 sage: from mjo.eja.eja_algebra import (HadamardEJA,
812 ....: RealSymmetricEJA,
813 ....: ComplexHermitianEJA,
814 ....: QuaternionHermitianEJA,
819 The rank of the Jordan spin algebra is always two::
821 sage: JordanSpinEJA(2).rank()
823 sage: JordanSpinEJA(3).rank()
825 sage: JordanSpinEJA(4).rank()
828 The rank of the `n`-by-`n` Hermitian real, complex, or
829 quaternion matrices is `n`::
831 sage: RealSymmetricEJA(4).rank()
833 sage: ComplexHermitianEJA(3).rank()
835 sage: QuaternionHermitianEJA(2).rank()
840 Ensure that every EJA that we know how to construct has a
841 positive integer rank, unless the algebra is trivial in
842 which case its rank will be zero::
844 sage: set_random_seed()
845 sage: J = random_eja()
849 sage: r > 0 or (r == 0 and J.is_trivial())
852 Ensure that computing the rank actually works, since the ranks
853 of all simple algebras are known and will be cached by default::
855 sage: J = HadamardEJA(4)
856 sage: J.rank.clear_cache()
862 sage: J = JordanSpinEJA(4)
863 sage: J.rank.clear_cache()
869 sage: J = RealSymmetricEJA(3)
870 sage: J.rank.clear_cache()
876 sage: J = ComplexHermitianEJA(2)
877 sage: J.rank.clear_cache()
883 sage: J = QuaternionHermitianEJA(2)
884 sage: J.rank.clear_cache()
895 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
896 R
= PolynomialRing(self
.base_ring(), var_names
)
900 # From a result in my book, these are the entries of the
901 # basis representation of L_x.
902 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
905 L_x
= matrix(R
, n
, n
, L_x_i_j
)
906 x_powers
= [ vars[k
]*(L_x
**k
)*self
.one().to_vector()
910 M
= matrix([x_powers
[0]])
914 M
= matrix(M
.rows() + [x_powers
[d
]])
916 # TODO: we've basically solved the system here.
917 # We should save the echelonized matrix somehow
918 # so that it can be reused in the charpoly method.
920 if new_rank
== old_rank
:
928 def vector_space(self
):
930 Return the vector space that underlies this algebra.
934 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
938 sage: J = RealSymmetricEJA(2)
939 sage: J.vector_space()
940 Vector space of dimension 3 over...
943 return self
.zero().to_vector().parent().ambient_vector_space()
946 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
949 class KnownRankEJA(object):
951 A class for algebras that we actually know we can construct. The
952 main issue is that, for most of our methods to make sense, we need
953 to know the rank of our algebra. Thus we can't simply generate a
954 "random" algebra, or even check that a given basis and product
955 satisfy the axioms; because even if everything looks OK, we wouldn't
956 know the rank we need to actuallty build the thing.
958 Not really a subclass of FDEJA because doing that causes method
959 resolution errors, e.g.
961 TypeError: Error when calling the metaclass bases
962 Cannot create a consistent method resolution
963 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
968 def _max_test_case_size():
970 Return an integer "size" that is an upper bound on the size of
971 this algebra when it is used in a random test
972 case. Unfortunately, the term "size" is quite vague -- when
973 dealing with `R^n` under either the Hadamard or Jordan spin
974 product, the "size" refers to the dimension `n`. When dealing
975 with a matrix algebra (real symmetric or complex/quaternion
976 Hermitian), it refers to the size of the matrix, which is
977 far less than the dimension of the underlying vector space.
979 We default to five in this class, which is safe in `R^n`. The
980 matrix algebra subclasses (or any class where the "size" is
981 interpreted to be far less than the dimension) should override
982 with a smaller number.
987 def random_instance(cls
, field
=AA
, **kwargs
):
989 Return a random instance of this type of algebra.
991 Beware, this will crash for "most instances" because the
992 constructor below looks wrong.
994 if cls
is TrivialEJA
:
995 # The TrivialEJA class doesn't take an "n" argument because
999 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
1000 return cls(n
, field
, **kwargs
)
1003 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1005 Return the Euclidean Jordan Algebra corresponding to the set
1006 `R^n` under the Hadamard product.
1008 Note: this is nothing more than the Cartesian product of ``n``
1009 copies of the spin algebra. Once Cartesian product algebras
1010 are implemented, this can go.
1014 sage: from mjo.eja.eja_algebra import HadamardEJA
1018 This multiplication table can be verified by hand::
1020 sage: J = HadamardEJA(3)
1021 sage: e0,e1,e2 = J.gens()
1037 We can change the generator prefix::
1039 sage: HadamardEJA(3, prefix='r').gens()
1043 def __init__(self
, n
, field
=AA
, **kwargs
):
1044 V
= VectorSpace(field
, n
)
1045 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
1048 fdeja
= super(HadamardEJA
, self
)
1049 fdeja
.__init
__(field
, mult_table
, **kwargs
)
1050 self
.rank
.set_cache(n
)
1052 def inner_product(self
, x
, y
):
1054 Faster to reimplement than to use natural representations.
1058 sage: from mjo.eja.eja_algebra import HadamardEJA
1062 Ensure that this is the usual inner product for the algebras
1065 sage: set_random_seed()
1066 sage: J = HadamardEJA.random_instance()
1067 sage: x,y = J.random_elements(2)
1068 sage: X = x.natural_representation()
1069 sage: Y = y.natural_representation()
1070 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1074 return x
.to_vector().inner_product(y
.to_vector())
1077 def random_eja(field
=AA
, nontrivial
=False):
1079 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1083 sage: from mjo.eja.eja_algebra import random_eja
1088 Euclidean Jordan algebra of dimension...
1091 eja_classes
= [HadamardEJA
,
1094 ComplexHermitianEJA
,
1095 QuaternionHermitianEJA
]
1097 eja_classes
.append(TrivialEJA
)
1098 classname
= choice(eja_classes
)
1099 return classname
.random_instance(field
=field
)
1106 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1108 def _max_test_case_size():
1109 # Play it safe, since this will be squared and the underlying
1110 # field can have dimension 4 (quaternions) too.
1113 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1115 Compared to the superclass constructor, we take a basis instead of
1116 a multiplication table because the latter can be computed in terms
1117 of the former when the product is known (like it is here).
1119 # Used in this class's fast _charpoly_coeff() override.
1120 self
._basis
_normalizers
= None
1122 # We're going to loop through this a few times, so now's a good
1123 # time to ensure that it isn't a generator expression.
1124 basis
= tuple(basis
)
1126 if len(basis
) > 1 and normalize_basis
:
1127 # We'll need sqrt(2) to normalize the basis, and this
1128 # winds up in the multiplication table, so the whole
1129 # algebra needs to be over the field extension.
1130 R
= PolynomialRing(field
, 'z')
1133 if p
.is_irreducible():
1134 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1135 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1136 self
._basis
_normalizers
= tuple(
1137 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1138 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1140 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1142 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
1143 fdeja
.__init
__(field
, Qs
, natural_basis
=basis
, **kwargs
)
1150 Override the parent method with something that tries to compute
1151 over a faster (non-extension) field.
1153 if self
._basis
_normalizers
is None:
1154 # We didn't normalize, so assume that the basis we started
1155 # with had entries in a nice field.
1156 return super(MatrixEuclideanJordanAlgebra
, self
).rank()
1158 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1159 self
._basis
_normalizers
) )
1161 # Do this over the rationals and convert back at the end.
1162 # Only works because we know the entries of the basis are
1164 J
= MatrixEuclideanJordanAlgebra(QQ
,
1166 normalize_basis
=False)
1170 def _charpoly_coeff(self
, i
):
1172 Override the parent method with something that tries to compute
1173 over a faster (non-extension) field.
1175 if self
._basis
_normalizers
is None:
1176 # We didn't normalize, so assume that the basis we started
1177 # with had entries in a nice field.
1178 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
1180 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1181 self
._basis
_normalizers
) )
1183 # Do this over the rationals and convert back at the end.
1184 J
= MatrixEuclideanJordanAlgebra(QQ
,
1186 normalize_basis
=False)
1187 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
1188 p
= J
._charpoly
_coeff
(i
)
1189 # p might be missing some vars, have to substitute "optionally"
1190 pairs
= zip(x
.base_ring().gens(), self
._basis
_normalizers
)
1191 substitutions
= { v: v*c for (v,c) in pairs }
1192 result
= p
.subs(substitutions
)
1194 # The result of "subs" can be either a coefficient-ring
1195 # element or a polynomial. Gotta handle both cases.
1197 return self
.base_ring()(result
)
1199 return result
.change_ring(self
.base_ring())
1203 def multiplication_table_from_matrix_basis(basis
):
1205 At least three of the five simple Euclidean Jordan algebras have the
1206 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1207 multiplication on the right is matrix multiplication. Given a basis
1208 for the underlying matrix space, this function returns a
1209 multiplication table (obtained by looping through the basis
1210 elements) for an algebra of those matrices.
1212 # In S^2, for example, we nominally have four coordinates even
1213 # though the space is of dimension three only. The vector space V
1214 # is supposed to hold the entire long vector, and the subspace W
1215 # of V will be spanned by the vectors that arise from symmetric
1216 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1217 field
= basis
[0].base_ring()
1218 dimension
= basis
[0].nrows()
1220 V
= VectorSpace(field
, dimension
**2)
1221 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1223 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1226 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1227 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1235 Embed the matrix ``M`` into a space of real matrices.
1237 The matrix ``M`` can have entries in any field at the moment:
1238 the real numbers, complex numbers, or quaternions. And although
1239 they are not a field, we can probably support octonions at some
1240 point, too. This function returns a real matrix that "acts like"
1241 the original with respect to matrix multiplication; i.e.
1243 real_embed(M*N) = real_embed(M)*real_embed(N)
1246 raise NotImplementedError
1250 def real_unembed(M
):
1252 The inverse of :meth:`real_embed`.
1254 raise NotImplementedError
1258 def natural_inner_product(cls
,X
,Y
):
1259 Xu
= cls
.real_unembed(X
)
1260 Yu
= cls
.real_unembed(Y
)
1261 tr
= (Xu
*Yu
).trace()
1264 # It's real already.
1267 # Otherwise, try the thing that works for complex numbers; and
1268 # if that doesn't work, the thing that works for quaternions.
1270 return tr
.vector()[0] # real part, imag part is index 1
1271 except AttributeError:
1272 # A quaternions doesn't have a vector() method, but does
1273 # have coefficient_tuple() method that returns the
1274 # coefficients of 1, i, j, and k -- in that order.
1275 return tr
.coefficient_tuple()[0]
1278 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1282 The identity function, for embedding real matrices into real
1288 def real_unembed(M
):
1290 The identity function, for unembedding real matrices from real
1296 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1298 The rank-n simple EJA consisting of real symmetric n-by-n
1299 matrices, the usual symmetric Jordan product, and the trace inner
1300 product. It has dimension `(n^2 + n)/2` over the reals.
1304 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1308 sage: J = RealSymmetricEJA(2)
1309 sage: e0, e1, e2 = J.gens()
1317 In theory, our "field" can be any subfield of the reals::
1319 sage: RealSymmetricEJA(2, RDF)
1320 Euclidean Jordan algebra of dimension 3 over Real Double Field
1321 sage: RealSymmetricEJA(2, RR)
1322 Euclidean Jordan algebra of dimension 3 over Real Field with
1323 53 bits of precision
1327 The dimension of this algebra is `(n^2 + n) / 2`::
1329 sage: set_random_seed()
1330 sage: n_max = RealSymmetricEJA._max_test_case_size()
1331 sage: n = ZZ.random_element(1, n_max)
1332 sage: J = RealSymmetricEJA(n)
1333 sage: J.dimension() == (n^2 + n)/2
1336 The Jordan multiplication is what we think it is::
1338 sage: set_random_seed()
1339 sage: J = RealSymmetricEJA.random_instance()
1340 sage: x,y = J.random_elements(2)
1341 sage: actual = (x*y).natural_representation()
1342 sage: X = x.natural_representation()
1343 sage: Y = y.natural_representation()
1344 sage: expected = (X*Y + Y*X)/2
1345 sage: actual == expected
1347 sage: J(expected) == x*y
1350 We can change the generator prefix::
1352 sage: RealSymmetricEJA(3, prefix='q').gens()
1353 (q0, q1, q2, q3, q4, q5)
1355 Our natural basis is normalized with respect to the natural inner
1356 product unless we specify otherwise::
1358 sage: set_random_seed()
1359 sage: J = RealSymmetricEJA.random_instance()
1360 sage: all( b.norm() == 1 for b in J.gens() )
1363 Since our natural basis is normalized with respect to the natural
1364 inner product, and since we know that this algebra is an EJA, any
1365 left-multiplication operator's matrix will be symmetric because
1366 natural->EJA basis representation is an isometry and within the EJA
1367 the operator is self-adjoint by the Jordan axiom::
1369 sage: set_random_seed()
1370 sage: x = RealSymmetricEJA.random_instance().random_element()
1371 sage: x.operator().matrix().is_symmetric()
1376 def _denormalized_basis(cls
, n
, field
):
1378 Return a basis for the space of real symmetric n-by-n matrices.
1382 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1386 sage: set_random_seed()
1387 sage: n = ZZ.random_element(1,5)
1388 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1389 sage: all( M.is_symmetric() for M in B)
1393 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1397 for j
in range(i
+1):
1398 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1402 Sij
= Eij
+ Eij
.transpose()
1408 def _max_test_case_size():
1409 return 4 # Dimension 10
1412 def __init__(self
, n
, field
=AA
, **kwargs
):
1413 basis
= self
._denormalized
_basis
(n
, field
)
1414 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, **kwargs
)
1415 self
.rank
.set_cache(n
)
1418 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1422 Embed the n-by-n complex matrix ``M`` into the space of real
1423 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1424 bi` to the block matrix ``[[a,b],[-b,a]]``.
1428 sage: from mjo.eja.eja_algebra import \
1429 ....: ComplexMatrixEuclideanJordanAlgebra
1433 sage: F = QuadraticField(-1, 'I')
1434 sage: x1 = F(4 - 2*i)
1435 sage: x2 = F(1 + 2*i)
1438 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1439 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1448 Embedding is a homomorphism (isomorphism, in fact)::
1450 sage: set_random_seed()
1451 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1452 sage: n = ZZ.random_element(n_max)
1453 sage: F = QuadraticField(-1, 'I')
1454 sage: X = random_matrix(F, n)
1455 sage: Y = random_matrix(F, n)
1456 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1457 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1458 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1465 raise ValueError("the matrix 'M' must be square")
1467 # We don't need any adjoined elements...
1468 field
= M
.base_ring().base_ring()
1472 a
= z
.list()[0] # real part, I guess
1473 b
= z
.list()[1] # imag part, I guess
1474 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1476 return matrix
.block(field
, n
, blocks
)
1480 def real_unembed(M
):
1482 The inverse of _embed_complex_matrix().
1486 sage: from mjo.eja.eja_algebra import \
1487 ....: ComplexMatrixEuclideanJordanAlgebra
1491 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1492 ....: [-2, 1, -4, 3],
1493 ....: [ 9, 10, 11, 12],
1494 ....: [-10, 9, -12, 11] ])
1495 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1497 [ 10*I + 9 12*I + 11]
1501 Unembedding is the inverse of embedding::
1503 sage: set_random_seed()
1504 sage: F = QuadraticField(-1, 'I')
1505 sage: M = random_matrix(F, 3)
1506 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1507 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1513 raise ValueError("the matrix 'M' must be square")
1514 if not n
.mod(2).is_zero():
1515 raise ValueError("the matrix 'M' must be a complex embedding")
1517 # If "M" was normalized, its base ring might have roots
1518 # adjoined and they can stick around after unembedding.
1519 field
= M
.base_ring()
1520 R
= PolynomialRing(field
, 'z')
1523 # Sage doesn't know how to embed AA into QQbar, i.e. how
1524 # to adjoin sqrt(-1) to AA.
1527 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1530 # Go top-left to bottom-right (reading order), converting every
1531 # 2-by-2 block we see to a single complex element.
1533 for k
in range(n
/2):
1534 for j
in range(n
/2):
1535 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1536 if submat
[0,0] != submat
[1,1]:
1537 raise ValueError('bad on-diagonal submatrix')
1538 if submat
[0,1] != -submat
[1,0]:
1539 raise ValueError('bad off-diagonal submatrix')
1540 z
= submat
[0,0] + submat
[0,1]*i
1543 return matrix(F
, n
/2, elements
)
1547 def natural_inner_product(cls
,X
,Y
):
1549 Compute a natural inner product in this algebra directly from
1554 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1558 This gives the same answer as the slow, default method implemented
1559 in :class:`MatrixEuclideanJordanAlgebra`::
1561 sage: set_random_seed()
1562 sage: J = ComplexHermitianEJA.random_instance()
1563 sage: x,y = J.random_elements(2)
1564 sage: Xe = x.natural_representation()
1565 sage: Ye = y.natural_representation()
1566 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1567 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1568 sage: expected = (X*Y).trace().real()
1569 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1570 sage: actual == expected
1574 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1577 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1579 The rank-n simple EJA consisting of complex Hermitian n-by-n
1580 matrices over the real numbers, the usual symmetric Jordan product,
1581 and the real-part-of-trace inner product. It has dimension `n^2` over
1586 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1590 In theory, our "field" can be any subfield of the reals::
1592 sage: ComplexHermitianEJA(2, RDF)
1593 Euclidean Jordan algebra of dimension 4 over Real Double Field
1594 sage: ComplexHermitianEJA(2, RR)
1595 Euclidean Jordan algebra of dimension 4 over Real Field with
1596 53 bits of precision
1600 The dimension of this algebra is `n^2`::
1602 sage: set_random_seed()
1603 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1604 sage: n = ZZ.random_element(1, n_max)
1605 sage: J = ComplexHermitianEJA(n)
1606 sage: J.dimension() == n^2
1609 The Jordan multiplication is what we think it is::
1611 sage: set_random_seed()
1612 sage: J = ComplexHermitianEJA.random_instance()
1613 sage: x,y = J.random_elements(2)
1614 sage: actual = (x*y).natural_representation()
1615 sage: X = x.natural_representation()
1616 sage: Y = y.natural_representation()
1617 sage: expected = (X*Y + Y*X)/2
1618 sage: actual == expected
1620 sage: J(expected) == x*y
1623 We can change the generator prefix::
1625 sage: ComplexHermitianEJA(2, prefix='z').gens()
1628 Our natural basis is normalized with respect to the natural inner
1629 product unless we specify otherwise::
1631 sage: set_random_seed()
1632 sage: J = ComplexHermitianEJA.random_instance()
1633 sage: all( b.norm() == 1 for b in J.gens() )
1636 Since our natural basis is normalized with respect to the natural
1637 inner product, and since we know that this algebra is an EJA, any
1638 left-multiplication operator's matrix will be symmetric because
1639 natural->EJA basis representation is an isometry and within the EJA
1640 the operator is self-adjoint by the Jordan axiom::
1642 sage: set_random_seed()
1643 sage: x = ComplexHermitianEJA.random_instance().random_element()
1644 sage: x.operator().matrix().is_symmetric()
1650 def _denormalized_basis(cls
, n
, field
):
1652 Returns a basis for the space of complex Hermitian n-by-n matrices.
1654 Why do we embed these? Basically, because all of numerical linear
1655 algebra assumes that you're working with vectors consisting of `n`
1656 entries from a field and scalars from the same field. There's no way
1657 to tell SageMath that (for example) the vectors contain complex
1658 numbers, while the scalar field is real.
1662 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1666 sage: set_random_seed()
1667 sage: n = ZZ.random_element(1,5)
1668 sage: field = QuadraticField(2, 'sqrt2')
1669 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1670 sage: all( M.is_symmetric() for M in B)
1674 R
= PolynomialRing(field
, 'z')
1676 F
= field
.extension(z
**2 + 1, 'I')
1679 # This is like the symmetric case, but we need to be careful:
1681 # * We want conjugate-symmetry, not just symmetry.
1682 # * The diagonal will (as a result) be real.
1686 for j
in range(i
+1):
1687 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1689 Sij
= cls
.real_embed(Eij
)
1692 # The second one has a minus because it's conjugated.
1693 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1695 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1698 # Since we embedded these, we can drop back to the "field" that we
1699 # started with instead of the complex extension "F".
1700 return ( s
.change_ring(field
) for s
in S
)
1703 def __init__(self
, n
, field
=AA
, **kwargs
):
1704 basis
= self
._denormalized
_basis
(n
,field
)
1705 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1706 self
.rank
.set_cache(n
)
1709 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1713 Embed the n-by-n quaternion matrix ``M`` into the space of real
1714 matrices of size 4n-by-4n by first sending each quaternion entry `z
1715 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1716 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1721 sage: from mjo.eja.eja_algebra import \
1722 ....: QuaternionMatrixEuclideanJordanAlgebra
1726 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1727 sage: i,j,k = Q.gens()
1728 sage: x = 1 + 2*i + 3*j + 4*k
1729 sage: M = matrix(Q, 1, [[x]])
1730 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1736 Embedding is a homomorphism (isomorphism, in fact)::
1738 sage: set_random_seed()
1739 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1740 sage: n = ZZ.random_element(n_max)
1741 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1742 sage: X = random_matrix(Q, n)
1743 sage: Y = random_matrix(Q, n)
1744 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1745 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1746 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1751 quaternions
= M
.base_ring()
1754 raise ValueError("the matrix 'M' must be square")
1756 F
= QuadraticField(-1, 'I')
1761 t
= z
.coefficient_tuple()
1766 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1767 [-c
+ d
*i
, a
- b
*i
]])
1768 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1769 blocks
.append(realM
)
1771 # We should have real entries by now, so use the realest field
1772 # we've got for the return value.
1773 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1778 def real_unembed(M
):
1780 The inverse of _embed_quaternion_matrix().
1784 sage: from mjo.eja.eja_algebra import \
1785 ....: QuaternionMatrixEuclideanJordanAlgebra
1789 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1790 ....: [-2, 1, -4, 3],
1791 ....: [-3, 4, 1, -2],
1792 ....: [-4, -3, 2, 1]])
1793 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1794 [1 + 2*i + 3*j + 4*k]
1798 Unembedding is the inverse of embedding::
1800 sage: set_random_seed()
1801 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1802 sage: M = random_matrix(Q, 3)
1803 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1804 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1810 raise ValueError("the matrix 'M' must be square")
1811 if not n
.mod(4).is_zero():
1812 raise ValueError("the matrix 'M' must be a quaternion embedding")
1814 # Use the base ring of the matrix to ensure that its entries can be
1815 # multiplied by elements of the quaternion algebra.
1816 field
= M
.base_ring()
1817 Q
= QuaternionAlgebra(field
,-1,-1)
1820 # Go top-left to bottom-right (reading order), converting every
1821 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1824 for l
in range(n
/4):
1825 for m
in range(n
/4):
1826 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1827 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1828 if submat
[0,0] != submat
[1,1].conjugate():
1829 raise ValueError('bad on-diagonal submatrix')
1830 if submat
[0,1] != -submat
[1,0].conjugate():
1831 raise ValueError('bad off-diagonal submatrix')
1832 z
= submat
[0,0].real()
1833 z
+= submat
[0,0].imag()*i
1834 z
+= submat
[0,1].real()*j
1835 z
+= submat
[0,1].imag()*k
1838 return matrix(Q
, n
/4, elements
)
1842 def natural_inner_product(cls
,X
,Y
):
1844 Compute a natural inner product in this algebra directly from
1849 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1853 This gives the same answer as the slow, default method implemented
1854 in :class:`MatrixEuclideanJordanAlgebra`::
1856 sage: set_random_seed()
1857 sage: J = QuaternionHermitianEJA.random_instance()
1858 sage: x,y = J.random_elements(2)
1859 sage: Xe = x.natural_representation()
1860 sage: Ye = y.natural_representation()
1861 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1862 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1863 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1864 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1865 sage: actual == expected
1869 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1872 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1875 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1876 matrices, the usual symmetric Jordan product, and the
1877 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1882 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1886 In theory, our "field" can be any subfield of the reals::
1888 sage: QuaternionHermitianEJA(2, RDF)
1889 Euclidean Jordan algebra of dimension 6 over Real Double Field
1890 sage: QuaternionHermitianEJA(2, RR)
1891 Euclidean Jordan algebra of dimension 6 over Real Field with
1892 53 bits of precision
1896 The dimension of this algebra is `2*n^2 - n`::
1898 sage: set_random_seed()
1899 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1900 sage: n = ZZ.random_element(1, n_max)
1901 sage: J = QuaternionHermitianEJA(n)
1902 sage: J.dimension() == 2*(n^2) - n
1905 The Jordan multiplication is what we think it is::
1907 sage: set_random_seed()
1908 sage: J = QuaternionHermitianEJA.random_instance()
1909 sage: x,y = J.random_elements(2)
1910 sage: actual = (x*y).natural_representation()
1911 sage: X = x.natural_representation()
1912 sage: Y = y.natural_representation()
1913 sage: expected = (X*Y + Y*X)/2
1914 sage: actual == expected
1916 sage: J(expected) == x*y
1919 We can change the generator prefix::
1921 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1922 (a0, a1, a2, a3, a4, a5)
1924 Our natural basis is normalized with respect to the natural inner
1925 product unless we specify otherwise::
1927 sage: set_random_seed()
1928 sage: J = QuaternionHermitianEJA.random_instance()
1929 sage: all( b.norm() == 1 for b in J.gens() )
1932 Since our natural basis is normalized with respect to the natural
1933 inner product, and since we know that this algebra is an EJA, any
1934 left-multiplication operator's matrix will be symmetric because
1935 natural->EJA basis representation is an isometry and within the EJA
1936 the operator is self-adjoint by the Jordan axiom::
1938 sage: set_random_seed()
1939 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1940 sage: x.operator().matrix().is_symmetric()
1945 def _denormalized_basis(cls
, n
, field
):
1947 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1949 Why do we embed these? Basically, because all of numerical
1950 linear algebra assumes that you're working with vectors consisting
1951 of `n` entries from a field and scalars from the same field. There's
1952 no way to tell SageMath that (for example) the vectors contain
1953 complex numbers, while the scalar field is real.
1957 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1961 sage: set_random_seed()
1962 sage: n = ZZ.random_element(1,5)
1963 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1964 sage: all( M.is_symmetric() for M in B )
1968 Q
= QuaternionAlgebra(QQ
,-1,-1)
1971 # This is like the symmetric case, but we need to be careful:
1973 # * We want conjugate-symmetry, not just symmetry.
1974 # * The diagonal will (as a result) be real.
1978 for j
in range(i
+1):
1979 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1981 Sij
= cls
.real_embed(Eij
)
1984 # The second, third, and fourth ones have a minus
1985 # because they're conjugated.
1986 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1988 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1990 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1992 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1995 # Since we embedded these, we can drop back to the "field" that we
1996 # started with instead of the quaternion algebra "Q".
1997 return ( s
.change_ring(field
) for s
in S
)
2000 def __init__(self
, n
, field
=AA
, **kwargs
):
2001 basis
= self
._denormalized
_basis
(n
,field
)
2002 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
2003 self
.rank
.set_cache(n
)
2006 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
2008 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2009 with the half-trace inner product and jordan product ``x*y =
2010 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
2011 symmetric positive-definite "bilinear form" matrix. It has
2012 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
2013 when ``B`` is the identity matrix of order ``n-1``.
2017 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2018 ....: JordanSpinEJA)
2022 When no bilinear form is specified, the identity matrix is used,
2023 and the resulting algebra is the Jordan spin algebra::
2025 sage: J0 = BilinearFormEJA(3)
2026 sage: J1 = JordanSpinEJA(3)
2027 sage: J0.multiplication_table() == J0.multiplication_table()
2032 We can create a zero-dimensional algebra::
2034 sage: J = BilinearFormEJA(0)
2038 We can check the multiplication condition given in the Jordan, von
2039 Neumann, and Wigner paper (and also discussed on my "On the
2040 symmetry..." paper). Note that this relies heavily on the standard
2041 choice of basis, as does anything utilizing the bilinear form matrix::
2043 sage: set_random_seed()
2044 sage: n = ZZ.random_element(5)
2045 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2046 sage: B = M.transpose()*M
2047 sage: J = BilinearFormEJA(n, B=B)
2048 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2049 sage: V = J.vector_space()
2050 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2051 ....: for ei in eis ]
2052 sage: actual = [ sis[i]*sis[j]
2053 ....: for i in range(n-1)
2054 ....: for j in range(n-1) ]
2055 sage: expected = [ J.one() if i == j else J.zero()
2056 ....: for i in range(n-1)
2057 ....: for j in range(n-1) ]
2058 sage: actual == expected
2061 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
2063 self
._B
= matrix
.identity(field
, max(0,n
-1))
2067 V
= VectorSpace(field
, n
)
2068 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2077 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
2078 zbar
= y0
*xbar
+ x0
*ybar
2079 z
= V([z0
] + zbar
.list())
2080 mult_table
[i
][j
] = z
2082 # The rank of this algebra is two, unless we're in a
2083 # one-dimensional ambient space (because the rank is bounded
2084 # by the ambient dimension).
2085 fdeja
= super(BilinearFormEJA
, self
)
2086 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2087 self
.rank
.set_cache(min(n
,2))
2089 def inner_product(self
, x
, y
):
2091 Half of the trace inner product.
2093 This is defined so that the special case of the Jordan spin
2094 algebra gets the usual inner product.
2098 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2102 Ensure that this is one-half of the trace inner-product when
2103 the algebra isn't just the reals (when ``n`` isn't one). This
2104 is in Faraut and Koranyi, and also my "On the symmetry..."
2107 sage: set_random_seed()
2108 sage: n = ZZ.random_element(2,5)
2109 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2110 sage: B = M.transpose()*M
2111 sage: J = BilinearFormEJA(n, B=B)
2112 sage: x = J.random_element()
2113 sage: y = J.random_element()
2114 sage: x.inner_product(y) == (x*y).trace()/2
2118 xvec
= x
.to_vector()
2120 yvec
= y
.to_vector()
2122 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
2125 class JordanSpinEJA(BilinearFormEJA
):
2127 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2128 with the usual inner product and jordan product ``x*y =
2129 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2134 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2138 This multiplication table can be verified by hand::
2140 sage: J = JordanSpinEJA(4)
2141 sage: e0,e1,e2,e3 = J.gens()
2157 We can change the generator prefix::
2159 sage: JordanSpinEJA(2, prefix='B').gens()
2164 Ensure that we have the usual inner product on `R^n`::
2166 sage: set_random_seed()
2167 sage: J = JordanSpinEJA.random_instance()
2168 sage: x,y = J.random_elements(2)
2169 sage: X = x.natural_representation()
2170 sage: Y = y.natural_representation()
2171 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2175 def __init__(self
, n
, field
=AA
, **kwargs
):
2176 # This is a special case of the BilinearFormEJA with the identity
2177 # matrix as its bilinear form.
2178 return super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2181 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
2183 The trivial Euclidean Jordan algebra consisting of only a zero element.
2187 sage: from mjo.eja.eja_algebra import TrivialEJA
2191 sage: J = TrivialEJA()
2198 sage: 7*J.one()*12*J.one()
2200 sage: J.one().inner_product(J.one())
2202 sage: J.one().norm()
2204 sage: J.one().subalgebra_generated_by()
2205 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2210 def __init__(self
, field
=AA
, **kwargs
):
2212 fdeja
= super(TrivialEJA
, self
)
2213 # The rank is zero using my definition, namely the dimension of the
2214 # largest subalgebra generated by any element.
2215 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2216 self
.rank
.set_cache(0)