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
65 sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
69 By definition, Jordan multiplication commutes::
71 sage: set_random_seed()
72 sage: J = random_eja()
73 sage: x,y = J.random_elements(2)
79 The ``field`` we're given must be real::
81 sage: JordanSpinEJA(2,QQbar)
82 Traceback (most recent call last):
84 ValueError: field is not real
88 if not field
.is_subring(RR
):
89 # Note: this does return true for the real algebraic
90 # field, and any quadratic field where we've specified
92 raise ValueError('field is not real')
95 self
._natural
_basis
= natural_basis
98 category
= MagmaticAlgebras(field
).FiniteDimensional()
99 category
= category
.WithBasis().Unital()
101 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
103 range(len(mult_table
)),
106 self
.print_options(bracket
='')
108 # The multiplication table we're given is necessarily in terms
109 # of vectors, because we don't have an algebra yet for
110 # anything to be an element of. However, it's faster in the
111 # long run to have the multiplication table be in terms of
112 # algebra elements. We do this after calling the superclass
113 # constructor so that from_vector() knows what to do.
114 self
._multiplication
_table
= [
115 list(map(lambda x
: self
.from_vector(x
), ls
))
120 def _element_constructor_(self
, elt
):
122 Construct an element of this algebra from its natural
125 This gets called only after the parent element _call_ method
126 fails to find a coercion for the argument.
130 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
132 ....: RealSymmetricEJA)
136 The identity in `S^n` is converted to the identity in the EJA::
138 sage: J = RealSymmetricEJA(3)
139 sage: I = matrix.identity(QQ,3)
140 sage: J(I) == J.one()
143 This skew-symmetric matrix can't be represented in the EJA::
145 sage: J = RealSymmetricEJA(3)
146 sage: A = matrix(QQ,3, lambda i,j: i-j)
148 Traceback (most recent call last):
150 ArithmeticError: vector is not in free module
154 Ensure that we can convert any element of the two non-matrix
155 simple algebras (whose natural representations are their usual
156 vector representations) back and forth faithfully::
158 sage: set_random_seed()
159 sage: J = HadamardEJA.random_instance()
160 sage: x = J.random_element()
161 sage: J(x.to_vector().column()) == x
163 sage: J = JordanSpinEJA.random_instance()
164 sage: x = J.random_element()
165 sage: J(x.to_vector().column()) == x
169 msg
= "not a naturally-represented algebra element"
171 # The superclass implementation of random_element()
172 # needs to be able to coerce "0" into the algebra.
174 elif elt
in self
.base_ring():
175 # Ensure that no base ring -> algebra coercion is performed
176 # by this method. There's some stupidity in sage that would
177 # otherwise propagate to this method; for example, sage thinks
178 # that the integer 3 belongs to the space of 2-by-2 matrices.
179 raise ValueError(msg
)
181 natural_basis
= self
.natural_basis()
182 basis_space
= natural_basis
[0].matrix_space()
183 if elt
not in basis_space
:
184 raise ValueError(msg
)
186 # Thanks for nothing! Matrix spaces aren't vector spaces in
187 # Sage, so we have to figure out its natural-basis coordinates
188 # ourselves. We use the basis space's ring instead of the
189 # element's ring because the basis space might be an algebraic
190 # closure whereas the base ring of the 3-by-3 identity matrix
191 # could be QQ instead of QQbar.
192 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
193 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
194 coords
= W
.coordinate_vector(_mat2vec(elt
))
195 return self
.from_vector(coords
)
200 Return a string representation of ``self``.
204 sage: from mjo.eja.eja_algebra import JordanSpinEJA
208 Ensure that it says what we think it says::
210 sage: JordanSpinEJA(2, field=AA)
211 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
212 sage: JordanSpinEJA(3, field=RDF)
213 Euclidean Jordan algebra of dimension 3 over Real Double Field
216 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
217 return fmt
.format(self
.dimension(), self
.base_ring())
219 def product_on_basis(self
, i
, j
):
220 return self
._multiplication
_table
[i
][j
]
222 def _a_regular_element(self
):
224 Guess a regular element. Needed to compute the basis for our
225 characteristic polynomial coefficients.
229 sage: from mjo.eja.eja_algebra import random_eja
233 Ensure that this hacky method succeeds for every algebra that we
234 know how to construct::
236 sage: set_random_seed()
237 sage: J = random_eja()
238 sage: J._a_regular_element().is_regular()
243 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
244 if not z
.is_regular():
245 raise ValueError("don't know a regular element")
250 def _charpoly_basis_space(self
):
252 Return the vector space spanned by the basis used in our
253 characteristic polynomial coefficients. This is used not only to
254 compute those coefficients, but also any time we need to
255 evaluate the coefficients (like when we compute the trace or
258 z
= self
._a
_regular
_element
()
259 # Don't use the parent vector space directly here in case this
260 # happens to be a subalgebra. In that case, we would be e.g.
261 # two-dimensional but span_of_basis() would expect three
263 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
264 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
265 V1
= V
.span_of_basis( basis
)
266 b
= (V1
.basis() + V1
.complement().basis())
267 return V
.span_of_basis(b
)
272 def _charpoly_coeff(self
, i
):
274 Return the coefficient polynomial "a_{i}" of this algebra's
275 general characteristic polynomial.
277 Having this be a separate cached method lets us compute and
278 store the trace/determinant (a_{r-1} and a_{0} respectively)
279 separate from the entire characteristic polynomial.
281 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
282 R
= A_of_x
.base_ring()
287 # Guaranteed by theory
290 # Danger: the in-place modification is done for performance
291 # reasons (reconstructing a matrix with huge polynomial
292 # entries is slow), but I don't know how cached_method works,
293 # so it's highly possible that we're modifying some global
294 # list variable by reference, here. In other words, you
295 # probably shouldn't call this method twice on the same
296 # algebra, at the same time, in two threads
297 Ai_orig
= A_of_x
.column(i
)
298 A_of_x
.set_column(i
,xr
)
299 numerator
= A_of_x
.det()
300 A_of_x
.set_column(i
,Ai_orig
)
302 # We're relying on the theory here to ensure that each a_i is
303 # indeed back in R, and the added negative signs are to make
304 # the whole charpoly expression sum to zero.
305 return R(-numerator
/detA
)
309 def _charpoly_matrix_system(self
):
311 Compute the matrix whose entries A_ij are polynomials in
312 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
313 corresponding to `x^r` and the determinent of the matrix A =
314 [A_ij]. In other words, all of the fixed (cachable) data needed
315 to compute the coefficients of the characteristic polynomial.
320 # Turn my vector space into a module so that "vectors" can
321 # have multivatiate polynomial entries.
322 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
323 R
= PolynomialRing(self
.base_ring(), names
)
325 # Using change_ring() on the parent's vector space doesn't work
326 # here because, in a subalgebra, that vector space has a basis
327 # and change_ring() tries to bring the basis along with it. And
328 # that doesn't work unless the new ring is a PID, which it usually
332 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
336 # And figure out the "left multiplication by x" matrix in
339 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
340 for i
in range(n
) ] # don't recompute these!
342 ek
= self
.monomial(k
).to_vector()
344 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
345 for i
in range(n
) ) )
346 Lx
= matrix
.column(R
, lmbx_cols
)
348 # Now we can compute powers of x "symbolically"
349 x_powers
= [self
.one().to_vector(), x
]
350 for d
in range(2, r
+1):
351 x_powers
.append( Lx
*(x_powers
[-1]) )
353 idmat
= matrix
.identity(R
, n
)
355 W
= self
._charpoly
_basis
_space
()
356 W
= W
.change_ring(R
.fraction_field())
358 # Starting with the standard coordinates x = (X1,X2,...,Xn)
359 # and then converting the entries to W-coordinates allows us
360 # to pass in the standard coordinates to the charpoly and get
361 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
364 # W.coordinates(x^2) eval'd at (standard z-coords)
368 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
370 # We want the middle equivalent thing in our matrix, but use
371 # the first equivalent thing instead so that we can pass in
372 # standard coordinates.
373 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
374 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
375 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
376 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
380 def characteristic_polynomial(self
):
382 Return a characteristic polynomial that works for all elements
385 The resulting polynomial has `n+1` variables, where `n` is the
386 dimension of this algebra. The first `n` variables correspond to
387 the coordinates of an algebra element: when evaluated at the
388 coordinates of an algebra element with respect to a certain
389 basis, the result is a univariate polynomial (in the one
390 remaining variable ``t``), namely the characteristic polynomial
395 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
399 The characteristic polynomial in the spin algebra is given in
400 Alizadeh, Example 11.11::
402 sage: J = JordanSpinEJA(3)
403 sage: p = J.characteristic_polynomial(); p
404 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
405 sage: xvec = J.one().to_vector()
409 By definition, the characteristic polynomial is a monic
410 degree-zero polynomial in a rank-zero algebra. Note that
411 Cayley-Hamilton is indeed satisfied since the polynomial
412 ``1`` evaluates to the identity element of the algebra on
415 sage: J = TrivialEJA()
416 sage: J.characteristic_polynomial()
423 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
424 a
= [ self
._charpoly
_coeff
(i
) for i
in range(r
+1) ]
426 # We go to a bit of trouble here to reorder the
427 # indeterminates, so that it's easier to evaluate the
428 # characteristic polynomial at x's coordinates and get back
429 # something in terms of t, which is what we want.
431 S
= PolynomialRing(self
.base_ring(),'t')
433 S
= PolynomialRing(S
, R
.variable_names())
436 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
439 def inner_product(self
, x
, y
):
441 The inner product associated with this Euclidean Jordan algebra.
443 Defaults to the trace inner product, but can be overridden by
444 subclasses if they are sure that the necessary properties are
449 sage: from mjo.eja.eja_algebra import random_eja
453 Our inner product is "associative," which means the following for
454 a symmetric bilinear form::
456 sage: set_random_seed()
457 sage: J = random_eja()
458 sage: x,y,z = J.random_elements(3)
459 sage: (x*y).inner_product(z) == y.inner_product(x*z)
463 X
= x
.natural_representation()
464 Y
= y
.natural_representation()
465 return self
.natural_inner_product(X
,Y
)
468 def is_trivial(self
):
470 Return whether or not this algebra is trivial.
472 A trivial algebra contains only the zero element.
476 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
481 sage: J = ComplexHermitianEJA(3)
487 sage: J = TrivialEJA()
492 return self
.dimension() == 0
495 def multiplication_table(self
):
497 Return a visual representation of this algebra's multiplication
498 table (on basis elements).
502 sage: from mjo.eja.eja_algebra import JordanSpinEJA
506 sage: J = JordanSpinEJA(4)
507 sage: J.multiplication_table()
508 +----++----+----+----+----+
509 | * || e0 | e1 | e2 | e3 |
510 +====++====+====+====+====+
511 | e0 || e0 | e1 | e2 | e3 |
512 +----++----+----+----+----+
513 | e1 || e1 | e0 | 0 | 0 |
514 +----++----+----+----+----+
515 | e2 || e2 | 0 | e0 | 0 |
516 +----++----+----+----+----+
517 | e3 || e3 | 0 | 0 | e0 |
518 +----++----+----+----+----+
521 M
= list(self
._multiplication
_table
) # copy
522 for i
in range(len(M
)):
523 # M had better be "square"
524 M
[i
] = [self
.monomial(i
)] + M
[i
]
525 M
= [["*"] + list(self
.gens())] + M
526 return table(M
, header_row
=True, header_column
=True, frame
=True)
529 def natural_basis(self
):
531 Return a more-natural representation of this algebra's basis.
533 Every finite-dimensional Euclidean Jordan Algebra is a direct
534 sum of five simple algebras, four of which comprise Hermitian
535 matrices. This method returns the original "natural" basis
536 for our underlying vector space. (Typically, the natural basis
537 is used to construct the multiplication table in the first place.)
539 Note that this will always return a matrix. The standard basis
540 in `R^n` will be returned as `n`-by-`1` column matrices.
544 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
545 ....: RealSymmetricEJA)
549 sage: J = RealSymmetricEJA(2)
551 Finite family {0: e0, 1: e1, 2: e2}
552 sage: J.natural_basis()
554 [1 0] [ 0 0.7071067811865475?] [0 0]
555 [0 0], [0.7071067811865475? 0], [0 1]
560 sage: J = JordanSpinEJA(2)
562 Finite family {0: e0, 1: e1}
563 sage: J.natural_basis()
570 if self
._natural
_basis
is None:
571 M
= self
.natural_basis_space()
572 return tuple( M(b
.to_vector()) for b
in self
.basis() )
574 return self
._natural
_basis
577 def natural_basis_space(self
):
579 Return the matrix space in which this algebra's natural basis
582 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
583 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
585 return self
._natural
_basis
[0].matrix_space()
589 def natural_inner_product(X
,Y
):
591 Compute the inner product of two naturally-represented elements.
593 For example in the real symmetric matrix EJA, this will compute
594 the trace inner-product of two n-by-n symmetric matrices. The
595 default should work for the real cartesian product EJA, the
596 Jordan spin EJA, and the real symmetric matrices. The others
597 will have to be overridden.
599 return (X
.conjugate_transpose()*Y
).trace()
605 Return the unit element of this algebra.
609 sage: from mjo.eja.eja_algebra import (HadamardEJA,
614 sage: J = HadamardEJA(5)
616 e0 + e1 + e2 + e3 + e4
620 The identity element acts like the identity::
622 sage: set_random_seed()
623 sage: J = random_eja()
624 sage: x = J.random_element()
625 sage: J.one()*x == x and x*J.one() == x
628 The matrix of the unit element's operator is the identity::
630 sage: set_random_seed()
631 sage: J = random_eja()
632 sage: actual = J.one().operator().matrix()
633 sage: expected = matrix.identity(J.base_ring(), J.dimension())
634 sage: actual == expected
638 # We can brute-force compute the matrices of the operators
639 # that correspond to the basis elements of this algebra.
640 # If some linear combination of those basis elements is the
641 # algebra identity, then the same linear combination of
642 # their matrices has to be the identity matrix.
644 # Of course, matrices aren't vectors in sage, so we have to
645 # appeal to the "long vectors" isometry.
646 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
648 # Now we use basis linear algebra to find the coefficients,
649 # of the matrices-as-vectors-linear-combination, which should
650 # work for the original algebra basis too.
651 A
= matrix
.column(self
.base_ring(), oper_vecs
)
653 # We used the isometry on the left-hand side already, but we
654 # still need to do it for the right-hand side. Recall that we
655 # wanted something that summed to the identity matrix.
656 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
658 # Now if there's an identity element in the algebra, this should work.
659 coeffs
= A
.solve_right(b
)
660 return self
.linear_combination(zip(self
.gens(), coeffs
))
663 def peirce_decomposition(self
, c
):
665 The Peirce decomposition of this algebra relative to the
668 In the future, this can be extended to a complete system of
669 orthogonal idempotents.
673 - ``c`` -- an idempotent of this algebra.
677 A triple (J0, J5, J1) containing two subalgebras and one subspace
680 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
681 corresponding to the eigenvalue zero.
683 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
684 corresponding to the eigenvalue one-half.
686 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
687 corresponding to the eigenvalue one.
689 These are the only possible eigenspaces for that operator, and this
690 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
691 orthogonal, and are subalgebras of this algebra with the appropriate
696 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
700 The canonical example comes from the symmetric matrices, which
701 decompose into diagonal and off-diagonal parts::
703 sage: J = RealSymmetricEJA(3)
704 sage: C = matrix(QQ, [ [1,0,0],
708 sage: J0,J5,J1 = J.peirce_decomposition(c)
710 Euclidean Jordan algebra of dimension 1...
712 Vector space of degree 6 and dimension 2...
714 Euclidean Jordan algebra of dimension 3...
718 Every algebra decomposes trivially with respect to its identity
721 sage: set_random_seed()
722 sage: J = random_eja()
723 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
724 sage: J0.dimension() == 0 and J5.dimension() == 0
726 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
729 The identity elements in the two subalgebras are the
730 projections onto their respective subspaces of the
731 superalgebra's identity element::
733 sage: set_random_seed()
734 sage: J = random_eja()
735 sage: x = J.random_element()
736 sage: if not J.is_trivial():
737 ....: while x.is_nilpotent():
738 ....: x = J.random_element()
739 sage: c = x.subalgebra_idempotent()
740 sage: J0,J5,J1 = J.peirce_decomposition(c)
741 sage: J1(c) == J1.one()
743 sage: J0(J.one() - c) == J0.one()
747 if not c
.is_idempotent():
748 raise ValueError("element is not idempotent: %s" % c
)
750 # Default these to what they should be if they turn out to be
751 # trivial, because eigenspaces_left() won't return eigenvalues
752 # corresponding to trivial spaces (e.g. it returns only the
753 # eigenspace corresponding to lambda=1 if you take the
754 # decomposition relative to the identity element).
755 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
756 J0
= trivial
# eigenvalue zero
757 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
758 J1
= trivial
# eigenvalue one
760 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
761 if eigval
== ~
(self
.base_ring()(2)):
764 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
765 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
771 raise ValueError("unexpected eigenvalue: %s" % eigval
)
776 def a_jordan_frame(self
):
778 Generate a Jordan frame for this algebra.
780 This implementation is based on the so-called "central
781 orthogonal idempotents" implemented for (semisimple) centers
782 of SageMath ``FiniteDimensionalAlgebrasWithBasis``. Since all
783 Euclidean Jordan algebas are commutative (and thus equal to
784 their own centers) and semisimple, the method should work more
785 or less as implemented, if it ever worked in the first place.
786 (I don't know the justification for the original implementation.
789 How it works: we loop through the algebras generators, looking
790 for their eigenspaces. If there's more than one eigenspace,
791 and if they result in more than one subalgebra, then we split
792 those subalgebras recursively until we get to subalgebras of
793 dimension one (whose idempotent is the unit element). Why does
794 some generator have to produce at least two subalgebras? I
795 dunno. But it seems to work.
797 Beware that Koecher defines the "center" of a Jordan algebra to
798 be something else, because the usual definition is stupid in a
799 (necessarily commutative) Jordan algebra.
803 sage: from mjo.eja.eja_algebra import (random_eja,
809 A Jordan frame for the trivial algebra has to be empty
810 (zero-length) since its rank is zero. More to the point, there
811 are no non-zero idempotents in the trivial EJA. This does not
812 cause any problems so long as we adopt the convention that the
813 empty sum is zero, since then the sole element of the trivial
814 EJA has an (empty) spectral decomposition::
816 sage: J = TrivialEJA()
817 sage: J.a_jordan_frame()
820 A one-dimensional algebra has rank one (equal to its dimension),
821 and only one primitive idempotent, namely the algebra's unit
824 sage: J = JordanSpinEJA(1)
825 sage: J.a_jordan_frame()
830 sage: J = random_eja()
831 sage: c = J.a_jordan_frame()
832 sage: all( x^2 == x for x in c )
835 sage: all( c[i]*c[j] == c[i]*(i==j) for i in range(r)
836 ....: for j in range(r) )
840 if self
.dimension() == 0:
842 if self
.dimension() == 1:
845 for g
in self
.gens():
846 eigenpairs
= g
.operator().matrix().right_eigenspaces()
847 if len(eigenpairs
) >= 2:
849 for eigval
, eigspace
in eigenpairs
:
850 # Make sub-EJAs from the matrix eigenspaces...
851 sb
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
853 # This will fail if e.g. the eigenspace basis
854 # contains two elements and their product
855 # isn't a linear combination of the two of
856 # them (i.e. the generated EJA isn't actually
858 s
= FiniteDimensionalEuclideanJordanSubalgebra(self
, sb
)
859 subalgebras
.append(s
)
860 except ArithmeticError as e
:
861 if str(e
) == "vector is not in free module":
862 # Ignore only the "not a sub-EJA" error
865 if len(subalgebras
) >= 2:
866 # apply this method recursively.
867 return tuple( c
.superalgebra_element()
868 for subalgebra
in subalgebras
869 for c
in subalgebra
.a_jordan_frame() )
871 # If we got here, the algebra didn't decompose, at least not when we looked at
872 # the eigenspaces corresponding only to basis elements of the algebra. The
873 # implementation I stole says that this should work because of Schur's Lemma,
874 # so I personally blame Schur's Lemma if it does not.
875 raise Exception("Schur's Lemma didn't work!")
878 def random_elements(self
, count
):
880 Return ``count`` random elements as a tuple.
884 sage: from mjo.eja.eja_algebra import JordanSpinEJA
888 sage: J = JordanSpinEJA(3)
889 sage: x,y,z = J.random_elements(3)
890 sage: all( [ x in J, y in J, z in J ])
892 sage: len( J.random_elements(10) ) == 10
896 return tuple( self
.random_element() for idx
in range(count
) )
901 Return the rank of this EJA.
905 The author knows of no algorithm to compute the rank of an EJA
906 where only the multiplication table is known. In lieu of one, we
907 require the rank to be specified when the algebra is created,
908 and simply pass along that number here.
912 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
913 ....: RealSymmetricEJA,
914 ....: ComplexHermitianEJA,
915 ....: QuaternionHermitianEJA,
920 The rank of the Jordan spin algebra is always two::
922 sage: JordanSpinEJA(2).rank()
924 sage: JordanSpinEJA(3).rank()
926 sage: JordanSpinEJA(4).rank()
929 The rank of the `n`-by-`n` Hermitian real, complex, or
930 quaternion matrices is `n`::
932 sage: RealSymmetricEJA(4).rank()
934 sage: ComplexHermitianEJA(3).rank()
936 sage: QuaternionHermitianEJA(2).rank()
941 Ensure that every EJA that we know how to construct has a
942 positive integer rank, unless the algebra is trivial in
943 which case its rank will be zero::
945 sage: set_random_seed()
946 sage: J = random_eja()
950 sage: r > 0 or (r == 0 and J.is_trivial())
957 def vector_space(self
):
959 Return the vector space that underlies this algebra.
963 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
967 sage: J = RealSymmetricEJA(2)
968 sage: J.vector_space()
969 Vector space of dimension 3 over...
972 return self
.zero().to_vector().parent().ambient_vector_space()
975 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
978 class KnownRankEJA(object):
980 A class for algebras that we actually know we can construct. The
981 main issue is that, for most of our methods to make sense, we need
982 to know the rank of our algebra. Thus we can't simply generate a
983 "random" algebra, or even check that a given basis and product
984 satisfy the axioms; because even if everything looks OK, we wouldn't
985 know the rank we need to actuallty build the thing.
987 Not really a subclass of FDEJA because doing that causes method
988 resolution errors, e.g.
990 TypeError: Error when calling the metaclass bases
991 Cannot create a consistent method resolution
992 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
997 def _max_test_case_size():
999 Return an integer "size" that is an upper bound on the size of
1000 this algebra when it is used in a random test
1001 case. Unfortunately, the term "size" is quite vague -- when
1002 dealing with `R^n` under either the Hadamard or Jordan spin
1003 product, the "size" refers to the dimension `n`. When dealing
1004 with a matrix algebra (real symmetric or complex/quaternion
1005 Hermitian), it refers to the size of the matrix, which is
1006 far less than the dimension of the underlying vector space.
1008 We default to five in this class, which is safe in `R^n`. The
1009 matrix algebra subclasses (or any class where the "size" is
1010 interpreted to be far less than the dimension) should override
1011 with a smaller number.
1016 def random_instance(cls
, field
=AA
, **kwargs
):
1018 Return a random instance of this type of algebra.
1020 Beware, this will crash for "most instances" because the
1021 constructor below looks wrong.
1023 if cls
is TrivialEJA
:
1024 # The TrivialEJA class doesn't take an "n" argument because
1028 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
1029 return cls(n
, field
, **kwargs
)
1032 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1034 Return the Euclidean Jordan Algebra corresponding to the set
1035 `R^n` under the Hadamard product.
1037 Note: this is nothing more than the Cartesian product of ``n``
1038 copies of the spin algebra. Once Cartesian product algebras
1039 are implemented, this can go.
1043 sage: from mjo.eja.eja_algebra import HadamardEJA
1047 This multiplication table can be verified by hand::
1049 sage: J = HadamardEJA(3)
1050 sage: e0,e1,e2 = J.gens()
1066 We can change the generator prefix::
1068 sage: HadamardEJA(3, prefix='r').gens()
1072 def __init__(self
, n
, field
=AA
, **kwargs
):
1073 V
= VectorSpace(field
, n
)
1074 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
1077 fdeja
= super(HadamardEJA
, self
)
1078 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
1080 def inner_product(self
, x
, y
):
1082 Faster to reimplement than to use natural representations.
1086 sage: from mjo.eja.eja_algebra import HadamardEJA
1090 Ensure that this is the usual inner product for the algebras
1093 sage: set_random_seed()
1094 sage: J = HadamardEJA.random_instance()
1095 sage: x,y = J.random_elements(2)
1096 sage: X = x.natural_representation()
1097 sage: Y = y.natural_representation()
1098 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1102 return x
.to_vector().inner_product(y
.to_vector())
1105 def random_eja(field
=AA
, nontrivial
=False):
1107 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1111 sage: from mjo.eja.eja_algebra import random_eja
1116 Euclidean Jordan algebra of dimension...
1119 eja_classes
= KnownRankEJA
.__subclasses
__()
1121 eja_classes
.remove(TrivialEJA
)
1122 classname
= choice(eja_classes
)
1123 return classname
.random_instance(field
=field
)
1130 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1132 def _max_test_case_size():
1133 # Play it safe, since this will be squared and the underlying
1134 # field can have dimension 4 (quaternions) too.
1137 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
1139 Compared to the superclass constructor, we take a basis instead of
1140 a multiplication table because the latter can be computed in terms
1141 of the former when the product is known (like it is here).
1143 # Used in this class's fast _charpoly_coeff() override.
1144 self
._basis
_normalizers
= None
1146 # We're going to loop through this a few times, so now's a good
1147 # time to ensure that it isn't a generator expression.
1148 basis
= tuple(basis
)
1150 if rank
> 1 and normalize_basis
:
1151 # We'll need sqrt(2) to normalize the basis, and this
1152 # winds up in the multiplication table, so the whole
1153 # algebra needs to be over the field extension.
1154 R
= PolynomialRing(field
, 'z')
1157 if p
.is_irreducible():
1158 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1159 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1160 self
._basis
_normalizers
= tuple(
1161 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1162 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1164 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1166 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
1167 return fdeja
.__init
__(field
,
1170 natural_basis
=basis
,
1175 def _charpoly_coeff(self
, i
):
1177 Override the parent method with something that tries to compute
1178 over a faster (non-extension) field.
1180 if self
._basis
_normalizers
is None:
1181 # We didn't normalize, so assume that the basis we started
1182 # with had entries in a nice field.
1183 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
1185 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1186 self
._basis
_normalizers
) )
1188 # Do this over the rationals and convert back at the end.
1189 J
= MatrixEuclideanJordanAlgebra(QQ
,
1192 normalize_basis
=False)
1193 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
1194 p
= J
._charpoly
_coeff
(i
)
1195 # p might be missing some vars, have to substitute "optionally"
1196 pairs
= zip(x
.base_ring().gens(), self
._basis
_normalizers
)
1197 substitutions
= { v: v*c for (v,c) in pairs }
1198 result
= p
.subs(substitutions
)
1200 # The result of "subs" can be either a coefficient-ring
1201 # element or a polynomial. Gotta handle both cases.
1203 return self
.base_ring()(result
)
1205 return result
.change_ring(self
.base_ring())
1209 def multiplication_table_from_matrix_basis(basis
):
1211 At least three of the five simple Euclidean Jordan algebras have the
1212 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1213 multiplication on the right is matrix multiplication. Given a basis
1214 for the underlying matrix space, this function returns a
1215 multiplication table (obtained by looping through the basis
1216 elements) for an algebra of those matrices.
1218 # In S^2, for example, we nominally have four coordinates even
1219 # though the space is of dimension three only. The vector space V
1220 # is supposed to hold the entire long vector, and the subspace W
1221 # of V will be spanned by the vectors that arise from symmetric
1222 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1223 field
= basis
[0].base_ring()
1224 dimension
= basis
[0].nrows()
1226 V
= VectorSpace(field
, dimension
**2)
1227 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1229 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1232 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1233 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1241 Embed the matrix ``M`` into a space of real matrices.
1243 The matrix ``M`` can have entries in any field at the moment:
1244 the real numbers, complex numbers, or quaternions. And although
1245 they are not a field, we can probably support octonions at some
1246 point, too. This function returns a real matrix that "acts like"
1247 the original with respect to matrix multiplication; i.e.
1249 real_embed(M*N) = real_embed(M)*real_embed(N)
1252 raise NotImplementedError
1256 def real_unembed(M
):
1258 The inverse of :meth:`real_embed`.
1260 raise NotImplementedError
1264 def natural_inner_product(cls
,X
,Y
):
1265 Xu
= cls
.real_unembed(X
)
1266 Yu
= cls
.real_unembed(Y
)
1267 tr
= (Xu
*Yu
).trace()
1270 # It's real already.
1273 # Otherwise, try the thing that works for complex numbers; and
1274 # if that doesn't work, the thing that works for quaternions.
1276 return tr
.vector()[0] # real part, imag part is index 1
1277 except AttributeError:
1278 # A quaternions doesn't have a vector() method, but does
1279 # have coefficient_tuple() method that returns the
1280 # coefficients of 1, i, j, and k -- in that order.
1281 return tr
.coefficient_tuple()[0]
1284 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1288 The identity function, for embedding real matrices into real
1294 def real_unembed(M
):
1296 The identity function, for unembedding real matrices from real
1302 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1304 The rank-n simple EJA consisting of real symmetric n-by-n
1305 matrices, the usual symmetric Jordan product, and the trace inner
1306 product. It has dimension `(n^2 + n)/2` over the reals.
1310 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1314 sage: J = RealSymmetricEJA(2)
1315 sage: e0, e1, e2 = J.gens()
1323 In theory, our "field" can be any subfield of the reals::
1325 sage: RealSymmetricEJA(2, RDF)
1326 Euclidean Jordan algebra of dimension 3 over Real Double Field
1327 sage: RealSymmetricEJA(2, RR)
1328 Euclidean Jordan algebra of dimension 3 over Real Field with
1329 53 bits of precision
1333 The dimension of this algebra is `(n^2 + n) / 2`::
1335 sage: set_random_seed()
1336 sage: n_max = RealSymmetricEJA._max_test_case_size()
1337 sage: n = ZZ.random_element(1, n_max)
1338 sage: J = RealSymmetricEJA(n)
1339 sage: J.dimension() == (n^2 + n)/2
1342 The Jordan multiplication is what we think it is::
1344 sage: set_random_seed()
1345 sage: J = RealSymmetricEJA.random_instance()
1346 sage: x,y = J.random_elements(2)
1347 sage: actual = (x*y).natural_representation()
1348 sage: X = x.natural_representation()
1349 sage: Y = y.natural_representation()
1350 sage: expected = (X*Y + Y*X)/2
1351 sage: actual == expected
1353 sage: J(expected) == x*y
1356 We can change the generator prefix::
1358 sage: RealSymmetricEJA(3, prefix='q').gens()
1359 (q0, q1, q2, q3, q4, q5)
1361 Our natural basis is normalized with respect to the natural inner
1362 product unless we specify otherwise::
1364 sage: set_random_seed()
1365 sage: J = RealSymmetricEJA.random_instance()
1366 sage: all( b.norm() == 1 for b in J.gens() )
1369 Since our natural basis is normalized with respect to the natural
1370 inner product, and since we know that this algebra is an EJA, any
1371 left-multiplication operator's matrix will be symmetric because
1372 natural->EJA basis representation is an isometry and within the EJA
1373 the operator is self-adjoint by the Jordan axiom::
1375 sage: set_random_seed()
1376 sage: x = RealSymmetricEJA.random_instance().random_element()
1377 sage: x.operator().matrix().is_symmetric()
1382 def _denormalized_basis(cls
, n
, field
):
1384 Return a basis for the space of real symmetric n-by-n matrices.
1388 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1392 sage: set_random_seed()
1393 sage: n = ZZ.random_element(1,5)
1394 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1395 sage: all( M.is_symmetric() for M in B)
1399 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1403 for j
in range(i
+1):
1404 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1408 Sij
= Eij
+ Eij
.transpose()
1414 def _max_test_case_size():
1415 return 4 # Dimension 10
1418 def __init__(self
, n
, field
=AA
, **kwargs
):
1419 basis
= self
._denormalized
_basis
(n
, field
)
1420 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1423 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1427 Embed the n-by-n complex matrix ``M`` into the space of real
1428 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1429 bi` to the block matrix ``[[a,b],[-b,a]]``.
1433 sage: from mjo.eja.eja_algebra import \
1434 ....: ComplexMatrixEuclideanJordanAlgebra
1438 sage: F = QuadraticField(-1, 'I')
1439 sage: x1 = F(4 - 2*i)
1440 sage: x2 = F(1 + 2*i)
1443 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1444 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1453 Embedding is a homomorphism (isomorphism, in fact)::
1455 sage: set_random_seed()
1456 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1457 sage: n = ZZ.random_element(n_max)
1458 sage: F = QuadraticField(-1, 'I')
1459 sage: X = random_matrix(F, n)
1460 sage: Y = random_matrix(F, n)
1461 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1462 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1463 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1470 raise ValueError("the matrix 'M' must be square")
1472 # We don't need any adjoined elements...
1473 field
= M
.base_ring().base_ring()
1477 a
= z
.list()[0] # real part, I guess
1478 b
= z
.list()[1] # imag part, I guess
1479 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1481 return matrix
.block(field
, n
, blocks
)
1485 def real_unembed(M
):
1487 The inverse of _embed_complex_matrix().
1491 sage: from mjo.eja.eja_algebra import \
1492 ....: ComplexMatrixEuclideanJordanAlgebra
1496 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1497 ....: [-2, 1, -4, 3],
1498 ....: [ 9, 10, 11, 12],
1499 ....: [-10, 9, -12, 11] ])
1500 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1502 [ 10*I + 9 12*I + 11]
1506 Unembedding is the inverse of embedding::
1508 sage: set_random_seed()
1509 sage: F = QuadraticField(-1, 'I')
1510 sage: M = random_matrix(F, 3)
1511 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1512 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1518 raise ValueError("the matrix 'M' must be square")
1519 if not n
.mod(2).is_zero():
1520 raise ValueError("the matrix 'M' must be a complex embedding")
1522 # If "M" was normalized, its base ring might have roots
1523 # adjoined and they can stick around after unembedding.
1524 field
= M
.base_ring()
1525 R
= PolynomialRing(field
, 'z')
1528 # Sage doesn't know how to embed AA into QQbar, i.e. how
1529 # to adjoin sqrt(-1) to AA.
1532 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1535 # Go top-left to bottom-right (reading order), converting every
1536 # 2-by-2 block we see to a single complex element.
1538 for k
in range(n
/2):
1539 for j
in range(n
/2):
1540 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1541 if submat
[0,0] != submat
[1,1]:
1542 raise ValueError('bad on-diagonal submatrix')
1543 if submat
[0,1] != -submat
[1,0]:
1544 raise ValueError('bad off-diagonal submatrix')
1545 z
= submat
[0,0] + submat
[0,1]*i
1548 return matrix(F
, n
/2, elements
)
1552 def natural_inner_product(cls
,X
,Y
):
1554 Compute a natural inner product in this algebra directly from
1559 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1563 This gives the same answer as the slow, default method implemented
1564 in :class:`MatrixEuclideanJordanAlgebra`::
1566 sage: set_random_seed()
1567 sage: J = ComplexHermitianEJA.random_instance()
1568 sage: x,y = J.random_elements(2)
1569 sage: Xe = x.natural_representation()
1570 sage: Ye = y.natural_representation()
1571 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1572 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1573 sage: expected = (X*Y).trace().real()
1574 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1575 sage: actual == expected
1579 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1582 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1584 The rank-n simple EJA consisting of complex Hermitian n-by-n
1585 matrices over the real numbers, the usual symmetric Jordan product,
1586 and the real-part-of-trace inner product. It has dimension `n^2` over
1591 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1595 In theory, our "field" can be any subfield of the reals::
1597 sage: ComplexHermitianEJA(2, RDF)
1598 Euclidean Jordan algebra of dimension 4 over Real Double Field
1599 sage: ComplexHermitianEJA(2, RR)
1600 Euclidean Jordan algebra of dimension 4 over Real Field with
1601 53 bits of precision
1605 The dimension of this algebra is `n^2`::
1607 sage: set_random_seed()
1608 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1609 sage: n = ZZ.random_element(1, n_max)
1610 sage: J = ComplexHermitianEJA(n)
1611 sage: J.dimension() == n^2
1614 The Jordan multiplication is what we think it is::
1616 sage: set_random_seed()
1617 sage: J = ComplexHermitianEJA.random_instance()
1618 sage: x,y = J.random_elements(2)
1619 sage: actual = (x*y).natural_representation()
1620 sage: X = x.natural_representation()
1621 sage: Y = y.natural_representation()
1622 sage: expected = (X*Y + Y*X)/2
1623 sage: actual == expected
1625 sage: J(expected) == x*y
1628 We can change the generator prefix::
1630 sage: ComplexHermitianEJA(2, prefix='z').gens()
1633 Our natural basis is normalized with respect to the natural inner
1634 product unless we specify otherwise::
1636 sage: set_random_seed()
1637 sage: J = ComplexHermitianEJA.random_instance()
1638 sage: all( b.norm() == 1 for b in J.gens() )
1641 Since our natural basis is normalized with respect to the natural
1642 inner product, and since we know that this algebra is an EJA, any
1643 left-multiplication operator's matrix will be symmetric because
1644 natural->EJA basis representation is an isometry and within the EJA
1645 the operator is self-adjoint by the Jordan axiom::
1647 sage: set_random_seed()
1648 sage: x = ComplexHermitianEJA.random_instance().random_element()
1649 sage: x.operator().matrix().is_symmetric()
1655 def _denormalized_basis(cls
, n
, field
):
1657 Returns a basis for the space of complex Hermitian n-by-n matrices.
1659 Why do we embed these? Basically, because all of numerical linear
1660 algebra assumes that you're working with vectors consisting of `n`
1661 entries from a field and scalars from the same field. There's no way
1662 to tell SageMath that (for example) the vectors contain complex
1663 numbers, while the scalar field is real.
1667 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1671 sage: set_random_seed()
1672 sage: n = ZZ.random_element(1,5)
1673 sage: field = QuadraticField(2, 'sqrt2')
1674 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1675 sage: all( M.is_symmetric() for M in B)
1679 R
= PolynomialRing(field
, 'z')
1681 F
= field
.extension(z
**2 + 1, 'I')
1684 # This is like the symmetric case, but we need to be careful:
1686 # * We want conjugate-symmetry, not just symmetry.
1687 # * The diagonal will (as a result) be real.
1691 for j
in range(i
+1):
1692 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1694 Sij
= cls
.real_embed(Eij
)
1697 # The second one has a minus because it's conjugated.
1698 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1700 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1703 # Since we embedded these, we can drop back to the "field" that we
1704 # started with instead of the complex extension "F".
1705 return ( s
.change_ring(field
) for s
in S
)
1708 def __init__(self
, n
, field
=AA
, **kwargs
):
1709 basis
= self
._denormalized
_basis
(n
,field
)
1710 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1713 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1717 Embed the n-by-n quaternion matrix ``M`` into the space of real
1718 matrices of size 4n-by-4n by first sending each quaternion entry `z
1719 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1720 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1725 sage: from mjo.eja.eja_algebra import \
1726 ....: QuaternionMatrixEuclideanJordanAlgebra
1730 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1731 sage: i,j,k = Q.gens()
1732 sage: x = 1 + 2*i + 3*j + 4*k
1733 sage: M = matrix(Q, 1, [[x]])
1734 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1740 Embedding is a homomorphism (isomorphism, in fact)::
1742 sage: set_random_seed()
1743 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1744 sage: n = ZZ.random_element(n_max)
1745 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1746 sage: X = random_matrix(Q, n)
1747 sage: Y = random_matrix(Q, n)
1748 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1749 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1750 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1755 quaternions
= M
.base_ring()
1758 raise ValueError("the matrix 'M' must be square")
1760 F
= QuadraticField(-1, 'I')
1765 t
= z
.coefficient_tuple()
1770 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1771 [-c
+ d
*i
, a
- b
*i
]])
1772 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1773 blocks
.append(realM
)
1775 # We should have real entries by now, so use the realest field
1776 # we've got for the return value.
1777 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1782 def real_unembed(M
):
1784 The inverse of _embed_quaternion_matrix().
1788 sage: from mjo.eja.eja_algebra import \
1789 ....: QuaternionMatrixEuclideanJordanAlgebra
1793 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1794 ....: [-2, 1, -4, 3],
1795 ....: [-3, 4, 1, -2],
1796 ....: [-4, -3, 2, 1]])
1797 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1798 [1 + 2*i + 3*j + 4*k]
1802 Unembedding is the inverse of embedding::
1804 sage: set_random_seed()
1805 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1806 sage: M = random_matrix(Q, 3)
1807 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1808 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1814 raise ValueError("the matrix 'M' must be square")
1815 if not n
.mod(4).is_zero():
1816 raise ValueError("the matrix 'M' must be a quaternion embedding")
1818 # Use the base ring of the matrix to ensure that its entries can be
1819 # multiplied by elements of the quaternion algebra.
1820 field
= M
.base_ring()
1821 Q
= QuaternionAlgebra(field
,-1,-1)
1824 # Go top-left to bottom-right (reading order), converting every
1825 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1828 for l
in range(n
/4):
1829 for m
in range(n
/4):
1830 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1831 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1832 if submat
[0,0] != submat
[1,1].conjugate():
1833 raise ValueError('bad on-diagonal submatrix')
1834 if submat
[0,1] != -submat
[1,0].conjugate():
1835 raise ValueError('bad off-diagonal submatrix')
1836 z
= submat
[0,0].real()
1837 z
+= submat
[0,0].imag()*i
1838 z
+= submat
[0,1].real()*j
1839 z
+= submat
[0,1].imag()*k
1842 return matrix(Q
, n
/4, elements
)
1846 def natural_inner_product(cls
,X
,Y
):
1848 Compute a natural inner product in this algebra directly from
1853 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1857 This gives the same answer as the slow, default method implemented
1858 in :class:`MatrixEuclideanJordanAlgebra`::
1860 sage: set_random_seed()
1861 sage: J = QuaternionHermitianEJA.random_instance()
1862 sage: x,y = J.random_elements(2)
1863 sage: Xe = x.natural_representation()
1864 sage: Ye = y.natural_representation()
1865 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1866 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1867 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1868 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1869 sage: actual == expected
1873 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1876 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1879 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1880 matrices, the usual symmetric Jordan product, and the
1881 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1886 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1890 In theory, our "field" can be any subfield of the reals::
1892 sage: QuaternionHermitianEJA(2, RDF)
1893 Euclidean Jordan algebra of dimension 6 over Real Double Field
1894 sage: QuaternionHermitianEJA(2, RR)
1895 Euclidean Jordan algebra of dimension 6 over Real Field with
1896 53 bits of precision
1900 The dimension of this algebra is `2*n^2 - n`::
1902 sage: set_random_seed()
1903 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1904 sage: n = ZZ.random_element(1, n_max)
1905 sage: J = QuaternionHermitianEJA(n)
1906 sage: J.dimension() == 2*(n^2) - n
1909 The Jordan multiplication is what we think it is::
1911 sage: set_random_seed()
1912 sage: J = QuaternionHermitianEJA.random_instance()
1913 sage: x,y = J.random_elements(2)
1914 sage: actual = (x*y).natural_representation()
1915 sage: X = x.natural_representation()
1916 sage: Y = y.natural_representation()
1917 sage: expected = (X*Y + Y*X)/2
1918 sage: actual == expected
1920 sage: J(expected) == x*y
1923 We can change the generator prefix::
1925 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1926 (a0, a1, a2, a3, a4, a5)
1928 Our natural basis is normalized with respect to the natural inner
1929 product unless we specify otherwise::
1931 sage: set_random_seed()
1932 sage: J = QuaternionHermitianEJA.random_instance()
1933 sage: all( b.norm() == 1 for b in J.gens() )
1936 Since our natural basis is normalized with respect to the natural
1937 inner product, and since we know that this algebra is an EJA, any
1938 left-multiplication operator's matrix will be symmetric because
1939 natural->EJA basis representation is an isometry and within the EJA
1940 the operator is self-adjoint by the Jordan axiom::
1942 sage: set_random_seed()
1943 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1944 sage: x.operator().matrix().is_symmetric()
1949 def _denormalized_basis(cls
, n
, field
):
1951 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1953 Why do we embed these? Basically, because all of numerical
1954 linear algebra assumes that you're working with vectors consisting
1955 of `n` entries from a field and scalars from the same field. There's
1956 no way to tell SageMath that (for example) the vectors contain
1957 complex numbers, while the scalar field is real.
1961 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1965 sage: set_random_seed()
1966 sage: n = ZZ.random_element(1,5)
1967 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1968 sage: all( M.is_symmetric() for M in B )
1972 Q
= QuaternionAlgebra(QQ
,-1,-1)
1975 # This is like the symmetric case, but we need to be careful:
1977 # * We want conjugate-symmetry, not just symmetry.
1978 # * The diagonal will (as a result) be real.
1982 for j
in range(i
+1):
1983 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1985 Sij
= cls
.real_embed(Eij
)
1988 # The second, third, and fourth ones have a minus
1989 # because they're conjugated.
1990 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1992 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1994 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1996 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1999 # Since we embedded these, we can drop back to the "field" that we
2000 # started with instead of the quaternion algebra "Q".
2001 return ( s
.change_ring(field
) for s
in S
)
2004 def __init__(self
, n
, field
=AA
, **kwargs
):
2005 basis
= self
._denormalized
_basis
(n
,field
)
2006 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
2009 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
2011 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2012 with the half-trace inner product and jordan product ``x*y =
2013 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
2014 symmetric positive-definite "bilinear form" matrix. It has
2015 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
2016 when ``B`` is the identity matrix of order ``n-1``.
2020 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2021 ....: JordanSpinEJA)
2025 When no bilinear form is specified, the identity matrix is used,
2026 and the resulting algebra is the Jordan spin algebra::
2028 sage: J0 = BilinearFormEJA(3)
2029 sage: J1 = JordanSpinEJA(3)
2030 sage: J0.multiplication_table() == J0.multiplication_table()
2035 We can create a zero-dimensional algebra::
2037 sage: J = BilinearFormEJA(0)
2041 We can check the multiplication condition given in the Jordan, von
2042 Neumann, and Wigner paper (and also discussed on my "On the
2043 symmetry..." paper). Note that this relies heavily on the standard
2044 choice of basis, as does anything utilizing the bilinear form matrix::
2046 sage: set_random_seed()
2047 sage: n = ZZ.random_element(5)
2048 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2049 sage: B = M.transpose()*M
2050 sage: J = BilinearFormEJA(n, B=B)
2051 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2052 sage: V = J.vector_space()
2053 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2054 ....: for ei in eis ]
2055 sage: actual = [ sis[i]*sis[j]
2056 ....: for i in range(n-1)
2057 ....: for j in range(n-1) ]
2058 sage: expected = [ J.one() if i == j else J.zero()
2059 ....: for i in range(n-1)
2060 ....: for j in range(n-1) ]
2061 sage: actual == expected
2064 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
2066 self
._B
= matrix
.identity(field
, max(0,n
-1))
2070 V
= VectorSpace(field
, n
)
2071 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2080 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
2081 zbar
= y0
*xbar
+ x0
*ybar
2082 z
= V([z0
] + zbar
.list())
2083 mult_table
[i
][j
] = z
2085 # The rank of this algebra is two, unless we're in a
2086 # one-dimensional ambient space (because the rank is bounded
2087 # by the ambient dimension).
2088 fdeja
= super(BilinearFormEJA
, self
)
2089 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
2091 def inner_product(self
, x
, y
):
2093 Half of the trace inner product.
2095 This is defined so that the special case of the Jordan spin
2096 algebra gets the usual inner product.
2100 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2104 Ensure that this is one-half of the trace inner-product when
2105 the algebra isn't just the reals (when ``n`` isn't one). This
2106 is in Faraut and Koranyi, and also my "On the symmetry..."
2109 sage: set_random_seed()
2110 sage: n = ZZ.random_element(2,5)
2111 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2112 sage: B = M.transpose()*M
2113 sage: J = BilinearFormEJA(n, B=B)
2114 sage: x = J.random_element()
2115 sage: y = J.random_element()
2116 sage: x.inner_product(y) == (x*y).trace()/2
2120 xvec
= x
.to_vector()
2122 yvec
= y
.to_vector()
2124 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
2127 class JordanSpinEJA(BilinearFormEJA
):
2129 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2130 with the usual inner product and jordan product ``x*y =
2131 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2136 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2140 This multiplication table can be verified by hand::
2142 sage: J = JordanSpinEJA(4)
2143 sage: e0,e1,e2,e3 = J.gens()
2159 We can change the generator prefix::
2161 sage: JordanSpinEJA(2, prefix='B').gens()
2166 Ensure that we have the usual inner product on `R^n`::
2168 sage: set_random_seed()
2169 sage: J = JordanSpinEJA.random_instance()
2170 sage: x,y = J.random_elements(2)
2171 sage: X = x.natural_representation()
2172 sage: Y = y.natural_representation()
2173 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2177 def __init__(self
, n
, field
=AA
, **kwargs
):
2178 # This is a special case of the BilinearFormEJA with the identity
2179 # matrix as its bilinear form.
2180 return super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2183 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
2185 The trivial Euclidean Jordan algebra consisting of only a zero element.
2189 sage: from mjo.eja.eja_algebra import TrivialEJA
2193 sage: J = TrivialEJA()
2200 sage: 7*J.one()*12*J.one()
2202 sage: J.one().inner_product(J.one())
2204 sage: J.one().norm()
2206 sage: J.one().subalgebra_generated_by()
2207 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2212 def __init__(self
, field
=AA
, **kwargs
):
2214 fdeja
= super(TrivialEJA
, self
)
2215 # The rank is zero using my definition, namely the dimension of the
2216 # largest subalgebra generated by any element.
2217 return fdeja
.__init
__(field
, mult_table
, rank
=0, **kwargs
)