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 (
65 ....: FiniteDimensionalEuclideanJordanAlgebra,
71 By definition, Jordan multiplication commutes::
73 sage: set_random_seed()
74 sage: J = random_eja()
75 sage: x,y = J.random_elements(2)
81 The ``field`` we're given must be real with ``check=True``::
83 sage: JordanSpinEJA(2,QQbar)
84 Traceback (most recent call last):
86 ValueError: field is not real
88 The multiplication table must be square with ``check=True``::
90 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
91 Traceback (most recent call last):
93 ValueError: multiplication table is not square
97 if not field
.is_subring(RR
):
98 # Note: this does return true for the real algebraic
99 # field, and any quadratic field where we've specified
101 raise ValueError('field is not real')
103 self
._natural
_basis
= natural_basis
106 category
= MagmaticAlgebras(field
).FiniteDimensional()
107 category
= category
.WithBasis().Unital()
109 # The multiplication table had better be square
112 if not all( len(l
) == n
for l
in mult_table
):
113 raise ValueError("multiplication table is not square")
115 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
120 self
.print_options(bracket
='')
122 # The multiplication table we're given is necessarily in terms
123 # of vectors, because we don't have an algebra yet for
124 # anything to be an element of. However, it's faster in the
125 # long run to have the multiplication table be in terms of
126 # algebra elements. We do this after calling the superclass
127 # constructor so that from_vector() knows what to do.
128 self
._multiplication
_table
= [
129 list(map(lambda x
: self
.from_vector(x
), ls
))
134 if not self
._is
_commutative
():
135 raise ValueError("algebra is not commutative")
136 if not self
._is
_jordanian
():
137 raise ValueError("Jordan identity does not hold")
138 if not self
._inner
_product
_is
_associative
():
139 raise ValueError("inner product is not associative")
141 def _element_constructor_(self
, elt
):
143 Construct an element of this algebra from its natural
146 This gets called only after the parent element _call_ method
147 fails to find a coercion for the argument.
151 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
153 ....: RealSymmetricEJA)
157 The identity in `S^n` is converted to the identity in the EJA::
159 sage: J = RealSymmetricEJA(3)
160 sage: I = matrix.identity(QQ,3)
161 sage: J(I) == J.one()
164 This skew-symmetric matrix can't be represented in the EJA::
166 sage: J = RealSymmetricEJA(3)
167 sage: A = matrix(QQ,3, lambda i,j: i-j)
169 Traceback (most recent call last):
171 ArithmeticError: vector is not in free module
175 Ensure that we can convert any element of the two non-matrix
176 simple algebras (whose natural representations are their usual
177 vector representations) back and forth faithfully::
179 sage: set_random_seed()
180 sage: J = HadamardEJA.random_instance()
181 sage: x = J.random_element()
182 sage: J(x.to_vector().column()) == x
184 sage: J = JordanSpinEJA.random_instance()
185 sage: x = J.random_element()
186 sage: J(x.to_vector().column()) == x
190 msg
= "not a naturally-represented algebra element"
192 # The superclass implementation of random_element()
193 # needs to be able to coerce "0" into the algebra.
195 elif elt
in self
.base_ring():
196 # Ensure that no base ring -> algebra coercion is performed
197 # by this method. There's some stupidity in sage that would
198 # otherwise propagate to this method; for example, sage thinks
199 # that the integer 3 belongs to the space of 2-by-2 matrices.
200 raise ValueError(msg
)
202 natural_basis
= self
.natural_basis()
203 basis_space
= natural_basis
[0].matrix_space()
204 if elt
not in basis_space
:
205 raise ValueError(msg
)
207 # Thanks for nothing! Matrix spaces aren't vector spaces in
208 # Sage, so we have to figure out its natural-basis coordinates
209 # ourselves. We use the basis space's ring instead of the
210 # element's ring because the basis space might be an algebraic
211 # closure whereas the base ring of the 3-by-3 identity matrix
212 # could be QQ instead of QQbar.
213 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
214 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
215 coords
= W
.coordinate_vector(_mat2vec(elt
))
216 return self
.from_vector(coords
)
219 def _max_test_case_size():
221 Return an integer "size" that is an upper bound on the size of
222 this algebra when it is used in a random test
223 case. Unfortunately, the term "size" is quite vague -- when
224 dealing with `R^n` under either the Hadamard or Jordan spin
225 product, the "size" refers to the dimension `n`. When dealing
226 with a matrix algebra (real symmetric or complex/quaternion
227 Hermitian), it refers to the size of the matrix, which is
228 far less than the dimension of the underlying vector space.
230 We default to five in this class, which is safe in `R^n`. The
231 matrix algebra subclasses (or any class where the "size" is
232 interpreted to be far less than the dimension) should override
233 with a smaller number.
239 Return a string representation of ``self``.
243 sage: from mjo.eja.eja_algebra import JordanSpinEJA
247 Ensure that it says what we think it says::
249 sage: JordanSpinEJA(2, field=AA)
250 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
251 sage: JordanSpinEJA(3, field=RDF)
252 Euclidean Jordan algebra of dimension 3 over Real Double Field
255 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
256 return fmt
.format(self
.dimension(), self
.base_ring())
258 def product_on_basis(self
, i
, j
):
259 return self
._multiplication
_table
[i
][j
]
261 def _is_commutative(self
):
263 Whether or not this algebra's multiplication table is commutative.
265 This method should of course always return ``True``, unless
266 this algebra was constructed with ``check=False`` and passed
267 an invalid multiplication table.
269 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
270 for i
in range(self
.dimension())
271 for j
in range(self
.dimension()) )
273 def _is_jordanian(self
):
275 Whether or not this algebra's multiplication table respects the
276 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
278 We only check one arrangement of `x` and `y`, so for a
279 ``True`` result to be truly true, you should also check
280 :meth:`_is_commutative`. This method should of course always
281 return ``True``, unless this algebra was constructed with
282 ``check=False`` and passed an invalid multiplication table.
284 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
286 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
287 for i
in range(self
.dimension())
288 for j
in range(self
.dimension()) )
290 def _inner_product_is_associative(self
):
292 Return whether or not this algebra's inner product `B` is
293 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
295 This method should of course always return ``True``, unless
296 this algebra was constructed with ``check=False`` and passed
297 an invalid multiplication table.
300 # Used to check whether or not something is zero in an inexact
301 # ring. This number is sufficient to allow the construction of
302 # QuaternionHermitianEJA(2, RDF) with check=True.
305 for i
in range(self
.dimension()):
306 for j
in range(self
.dimension()):
307 for k
in range(self
.dimension()):
311 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
313 if self
.base_ring().is_exact():
317 if diff
.abs() > epsilon
:
323 def characteristic_polynomial_of(self
):
325 Return the algebra's "characteristic polynomial of" function,
326 which is itself a multivariate polynomial that, when evaluated
327 at the coordinates of some algebra element, returns that
328 element's characteristic polynomial.
330 The resulting polynomial has `n+1` variables, where `n` is the
331 dimension of this algebra. The first `n` variables correspond to
332 the coordinates of an algebra element: when evaluated at the
333 coordinates of an algebra element with respect to a certain
334 basis, the result is a univariate polynomial (in the one
335 remaining variable ``t``), namely the characteristic polynomial
340 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
344 The characteristic polynomial in the spin algebra is given in
345 Alizadeh, Example 11.11::
347 sage: J = JordanSpinEJA(3)
348 sage: p = J.characteristic_polynomial_of(); p
349 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
350 sage: xvec = J.one().to_vector()
354 By definition, the characteristic polynomial is a monic
355 degree-zero polynomial in a rank-zero algebra. Note that
356 Cayley-Hamilton is indeed satisfied since the polynomial
357 ``1`` evaluates to the identity element of the algebra on
360 sage: J = TrivialEJA()
361 sage: J.characteristic_polynomial_of()
368 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
369 a
= self
._charpoly
_coefficients
()
371 # We go to a bit of trouble here to reorder the
372 # indeterminates, so that it's easier to evaluate the
373 # characteristic polynomial at x's coordinates and get back
374 # something in terms of t, which is what we want.
375 S
= PolynomialRing(self
.base_ring(),'t')
379 S
= PolynomialRing(S
, R
.variable_names())
382 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
385 def inner_product(self
, x
, y
):
387 The inner product associated with this Euclidean Jordan algebra.
389 Defaults to the trace inner product, but can be overridden by
390 subclasses if they are sure that the necessary properties are
395 sage: from mjo.eja.eja_algebra import random_eja
399 Our inner product is "associative," which means the following for
400 a symmetric bilinear form::
402 sage: set_random_seed()
403 sage: J = random_eja()
404 sage: x,y,z = J.random_elements(3)
405 sage: (x*y).inner_product(z) == y.inner_product(x*z)
409 X
= x
.natural_representation()
410 Y
= y
.natural_representation()
411 return self
.natural_inner_product(X
,Y
)
414 def is_trivial(self
):
416 Return whether or not this algebra is trivial.
418 A trivial algebra contains only the zero element.
422 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
427 sage: J = ComplexHermitianEJA(3)
433 sage: J = TrivialEJA()
438 return self
.dimension() == 0
441 def multiplication_table(self
):
443 Return a visual representation of this algebra's multiplication
444 table (on basis elements).
448 sage: from mjo.eja.eja_algebra import JordanSpinEJA
452 sage: J = JordanSpinEJA(4)
453 sage: J.multiplication_table()
454 +----++----+----+----+----+
455 | * || e0 | e1 | e2 | e3 |
456 +====++====+====+====+====+
457 | e0 || e0 | e1 | e2 | e3 |
458 +----++----+----+----+----+
459 | e1 || e1 | e0 | 0 | 0 |
460 +----++----+----+----+----+
461 | e2 || e2 | 0 | e0 | 0 |
462 +----++----+----+----+----+
463 | e3 || e3 | 0 | 0 | e0 |
464 +----++----+----+----+----+
467 M
= list(self
._multiplication
_table
) # copy
468 for i
in range(len(M
)):
469 # M had better be "square"
470 M
[i
] = [self
.monomial(i
)] + M
[i
]
471 M
= [["*"] + list(self
.gens())] + M
472 return table(M
, header_row
=True, header_column
=True, frame
=True)
475 def natural_basis(self
):
477 Return a more-natural representation of this algebra's basis.
479 Every finite-dimensional Euclidean Jordan Algebra is a direct
480 sum of five simple algebras, four of which comprise Hermitian
481 matrices. This method returns the original "natural" basis
482 for our underlying vector space. (Typically, the natural basis
483 is used to construct the multiplication table in the first place.)
485 Note that this will always return a matrix. The standard basis
486 in `R^n` will be returned as `n`-by-`1` column matrices.
490 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
491 ....: RealSymmetricEJA)
495 sage: J = RealSymmetricEJA(2)
497 Finite family {0: e0, 1: e1, 2: e2}
498 sage: J.natural_basis()
500 [1 0] [ 0 0.7071067811865475?] [0 0]
501 [0 0], [0.7071067811865475? 0], [0 1]
506 sage: J = JordanSpinEJA(2)
508 Finite family {0: e0, 1: e1}
509 sage: J.natural_basis()
516 if self
._natural
_basis
is None:
517 M
= self
.natural_basis_space()
518 return tuple( M(b
.to_vector()) for b
in self
.basis() )
520 return self
._natural
_basis
523 def natural_basis_space(self
):
525 Return the matrix space in which this algebra's natural basis
528 Generally this will be an `n`-by-`1` column-vector space,
529 except when the algebra is trivial. There it's `n`-by-`n`
530 (where `n` is zero), to ensure that two elements of the
531 natural basis space (empty matrices) can be multiplied.
533 if self
.is_trivial():
534 return MatrixSpace(self
.base_ring(), 0)
535 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
536 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
538 return self
._natural
_basis
[0].matrix_space()
542 def natural_inner_product(X
,Y
):
544 Compute the inner product of two naturally-represented elements.
546 For example in the real symmetric matrix EJA, this will compute
547 the trace inner-product of two n-by-n symmetric matrices. The
548 default should work for the real cartesian product EJA, the
549 Jordan spin EJA, and the real symmetric matrices. The others
550 will have to be overridden.
552 return (X
.conjugate_transpose()*Y
).trace()
558 Return the unit element of this algebra.
562 sage: from mjo.eja.eja_algebra import (HadamardEJA,
567 sage: J = HadamardEJA(5)
569 e0 + e1 + e2 + e3 + e4
573 The identity element acts like the identity::
575 sage: set_random_seed()
576 sage: J = random_eja()
577 sage: x = J.random_element()
578 sage: J.one()*x == x and x*J.one() == x
581 The matrix of the unit element's operator is the identity::
583 sage: set_random_seed()
584 sage: J = random_eja()
585 sage: actual = J.one().operator().matrix()
586 sage: expected = matrix.identity(J.base_ring(), J.dimension())
587 sage: actual == expected
591 # We can brute-force compute the matrices of the operators
592 # that correspond to the basis elements of this algebra.
593 # If some linear combination of those basis elements is the
594 # algebra identity, then the same linear combination of
595 # their matrices has to be the identity matrix.
597 # Of course, matrices aren't vectors in sage, so we have to
598 # appeal to the "long vectors" isometry.
599 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
601 # Now we use basis linear algebra to find the coefficients,
602 # of the matrices-as-vectors-linear-combination, which should
603 # work for the original algebra basis too.
604 A
= matrix
.column(self
.base_ring(), oper_vecs
)
606 # We used the isometry on the left-hand side already, but we
607 # still need to do it for the right-hand side. Recall that we
608 # wanted something that summed to the identity matrix.
609 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
611 # Now if there's an identity element in the algebra, this should work.
612 coeffs
= A
.solve_right(b
)
613 return self
.linear_combination(zip(self
.gens(), coeffs
))
616 def peirce_decomposition(self
, c
):
618 The Peirce decomposition of this algebra relative to the
621 In the future, this can be extended to a complete system of
622 orthogonal idempotents.
626 - ``c`` -- an idempotent of this algebra.
630 A triple (J0, J5, J1) containing two subalgebras and one subspace
633 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
634 corresponding to the eigenvalue zero.
636 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
637 corresponding to the eigenvalue one-half.
639 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
640 corresponding to the eigenvalue one.
642 These are the only possible eigenspaces for that operator, and this
643 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
644 orthogonal, and are subalgebras of this algebra with the appropriate
649 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
653 The canonical example comes from the symmetric matrices, which
654 decompose into diagonal and off-diagonal parts::
656 sage: J = RealSymmetricEJA(3)
657 sage: C = matrix(QQ, [ [1,0,0],
661 sage: J0,J5,J1 = J.peirce_decomposition(c)
663 Euclidean Jordan algebra of dimension 1...
665 Vector space of degree 6 and dimension 2...
667 Euclidean Jordan algebra of dimension 3...
668 sage: J0.one().natural_representation()
672 sage: orig_df = AA.options.display_format
673 sage: AA.options.display_format = 'radical'
674 sage: J.from_vector(J5.basis()[0]).natural_representation()
678 sage: J.from_vector(J5.basis()[1]).natural_representation()
682 sage: AA.options.display_format = orig_df
683 sage: J1.one().natural_representation()
690 Every algebra decomposes trivially with respect to its identity
693 sage: set_random_seed()
694 sage: J = random_eja()
695 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
696 sage: J0.dimension() == 0 and J5.dimension() == 0
698 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
701 The decomposition is into eigenspaces, and its components are
702 therefore necessarily orthogonal. Moreover, the identity
703 elements in the two subalgebras are the projections onto their
704 respective subspaces of the superalgebra's identity element::
706 sage: set_random_seed()
707 sage: J = random_eja()
708 sage: x = J.random_element()
709 sage: if not J.is_trivial():
710 ....: while x.is_nilpotent():
711 ....: x = J.random_element()
712 sage: c = x.subalgebra_idempotent()
713 sage: J0,J5,J1 = J.peirce_decomposition(c)
715 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
716 ....: w = w.superalgebra_element()
717 ....: y = J.from_vector(y)
718 ....: z = z.superalgebra_element()
719 ....: ipsum += w.inner_product(y).abs()
720 ....: ipsum += w.inner_product(z).abs()
721 ....: ipsum += y.inner_product(z).abs()
724 sage: J1(c) == J1.one()
726 sage: J0(J.one() - c) == J0.one()
730 if not c
.is_idempotent():
731 raise ValueError("element is not idempotent: %s" % c
)
733 # Default these to what they should be if they turn out to be
734 # trivial, because eigenspaces_left() won't return eigenvalues
735 # corresponding to trivial spaces (e.g. it returns only the
736 # eigenspace corresponding to lambda=1 if you take the
737 # decomposition relative to the identity element).
738 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
739 J0
= trivial
# eigenvalue zero
740 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
741 J1
= trivial
# eigenvalue one
743 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
744 if eigval
== ~
(self
.base_ring()(2)):
747 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
748 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
754 raise ValueError("unexpected eigenvalue: %s" % eigval
)
759 def random_element(self
, thorough
=False):
761 Return a random element of this algebra.
763 Our algebra superclass method only returns a linear
764 combination of at most two basis elements. We instead
765 want the vector space "random element" method that
766 returns a more diverse selection.
770 - ``thorough`` -- (boolean; default False) whether or not we
771 should generate irrational coefficients for the random
772 element when our base ring is irrational; this slows the
773 algebra operations to a crawl, but any truly random method
777 # For a general base ring... maybe we can trust this to do the
778 # right thing? Unlikely, but.
779 V
= self
.vector_space()
780 v
= V
.random_element()
782 if self
.base_ring() is AA
:
783 # The "random element" method of the algebraic reals is
784 # stupid at the moment, and only returns integers between
785 # -2 and 2, inclusive:
787 # https://trac.sagemath.org/ticket/30875
789 # Instead, we implement our own "random vector" method,
790 # and then coerce that into the algebra. We use the vector
791 # space degree here instead of the dimension because a
792 # subalgebra could (for example) be spanned by only two
793 # vectors, each with five coordinates. We need to
794 # generate all five coordinates.
796 v
*= QQbar
.random_element().real()
798 v
*= QQ
.random_element()
800 return self
.from_vector(V
.coordinate_vector(v
))
802 def random_elements(self
, count
, thorough
=False):
804 Return ``count`` random elements as a tuple.
808 - ``thorough`` -- (boolean; default False) whether or not we
809 should generate irrational coefficients for the random
810 elements when our base ring is irrational; this slows the
811 algebra operations to a crawl, but any truly random method
816 sage: from mjo.eja.eja_algebra import JordanSpinEJA
820 sage: J = JordanSpinEJA(3)
821 sage: x,y,z = J.random_elements(3)
822 sage: all( [ x in J, y in J, z in J ])
824 sage: len( J.random_elements(10) ) == 10
828 return tuple( self
.random_element(thorough
)
829 for idx
in range(count
) )
832 def random_instance(cls
, field
=AA
, **kwargs
):
834 Return a random instance of this type of algebra.
836 Beware, this will crash for "most instances" because the
837 constructor below looks wrong.
839 if cls
is TrivialEJA
:
840 # The TrivialEJA class doesn't take an "n" argument because
844 n
= ZZ
.random_element(cls
._max
_test
_case
_size
() + 1)
845 return cls(n
, field
, **kwargs
)
848 def _charpoly_coefficients(self
):
850 The `r` polynomial coefficients of the "characteristic polynomial
854 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
855 R
= PolynomialRing(self
.base_ring(), var_names
)
857 F
= R
.fraction_field()
860 # From a result in my book, these are the entries of the
861 # basis representation of L_x.
862 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
865 L_x
= matrix(F
, n
, n
, L_x_i_j
)
868 if self
.rank
.is_in_cache():
870 # There's no need to pad the system with redundant
871 # columns if we *know* they'll be redundant.
874 # Compute an extra power in case the rank is equal to
875 # the dimension (otherwise, we would stop at x^(r-1)).
876 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
877 for k
in range(n
+1) ]
878 A
= matrix
.column(F
, x_powers
[:n
])
879 AE
= A
.extended_echelon_form()
886 # The theory says that only the first "r" coefficients are
887 # nonzero, and they actually live in the original polynomial
888 # ring and not the fraction field. We negate them because
889 # in the actual characteristic polynomial, they get moved
890 # to the other side where x^r lives.
891 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
896 Return the rank of this EJA.
898 This is a cached method because we know the rank a priori for
899 all of the algebras we can construct. Thus we can avoid the
900 expensive ``_charpoly_coefficients()`` call unless we truly
901 need to compute the whole characteristic polynomial.
905 sage: from mjo.eja.eja_algebra import (HadamardEJA,
907 ....: RealSymmetricEJA,
908 ....: ComplexHermitianEJA,
909 ....: QuaternionHermitianEJA,
914 The rank of the Jordan spin algebra is always two::
916 sage: JordanSpinEJA(2).rank()
918 sage: JordanSpinEJA(3).rank()
920 sage: JordanSpinEJA(4).rank()
923 The rank of the `n`-by-`n` Hermitian real, complex, or
924 quaternion matrices is `n`::
926 sage: RealSymmetricEJA(4).rank()
928 sage: ComplexHermitianEJA(3).rank()
930 sage: QuaternionHermitianEJA(2).rank()
935 Ensure that every EJA that we know how to construct has a
936 positive integer rank, unless the algebra is trivial in
937 which case its rank will be zero::
939 sage: set_random_seed()
940 sage: J = random_eja()
944 sage: r > 0 or (r == 0 and J.is_trivial())
947 Ensure that computing the rank actually works, since the ranks
948 of all simple algebras are known and will be cached by default::
950 sage: J = HadamardEJA(4)
951 sage: J.rank.clear_cache()
957 sage: J = JordanSpinEJA(4)
958 sage: J.rank.clear_cache()
964 sage: J = RealSymmetricEJA(3)
965 sage: J.rank.clear_cache()
971 sage: J = ComplexHermitianEJA(2)
972 sage: J.rank.clear_cache()
978 sage: J = QuaternionHermitianEJA(2)
979 sage: J.rank.clear_cache()
983 return len(self
._charpoly
_coefficients
())
986 def vector_space(self
):
988 Return the vector space that underlies this algebra.
992 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
996 sage: J = RealSymmetricEJA(2)
997 sage: J.vector_space()
998 Vector space of dimension 3 over...
1001 return self
.zero().to_vector().parent().ambient_vector_space()
1004 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1007 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1009 Return the Euclidean Jordan Algebra corresponding to the set
1010 `R^n` under the Hadamard product.
1012 Note: this is nothing more than the Cartesian product of ``n``
1013 copies of the spin algebra. Once Cartesian product algebras
1014 are implemented, this can go.
1018 sage: from mjo.eja.eja_algebra import HadamardEJA
1022 This multiplication table can be verified by hand::
1024 sage: J = HadamardEJA(3)
1025 sage: e0,e1,e2 = J.gens()
1041 We can change the generator prefix::
1043 sage: HadamardEJA(3, prefix='r').gens()
1047 def __init__(self
, n
, field
=AA
, **kwargs
):
1048 V
= VectorSpace(field
, n
)
1049 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
1052 fdeja
= super(HadamardEJA
, self
)
1053 fdeja
.__init
__(field
, mult_table
, **kwargs
)
1054 self
.rank
.set_cache(n
)
1056 def inner_product(self
, x
, y
):
1058 Faster to reimplement than to use natural representations.
1062 sage: from mjo.eja.eja_algebra import HadamardEJA
1066 Ensure that this is the usual inner product for the algebras
1069 sage: set_random_seed()
1070 sage: J = HadamardEJA.random_instance()
1071 sage: x,y = J.random_elements(2)
1072 sage: X = x.natural_representation()
1073 sage: Y = y.natural_representation()
1074 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1078 return x
.to_vector().inner_product(y
.to_vector())
1081 def random_eja(field
=AA
):
1083 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1087 sage: from mjo.eja.eja_algebra import random_eja
1092 Euclidean Jordan algebra of dimension...
1095 classname
= choice([TrivialEJA
,
1099 ComplexHermitianEJA
,
1100 QuaternionHermitianEJA
])
1101 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_coefficients() 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
)
1148 def _charpoly_coefficients(self
):
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
)._charpoly
_coefficients
()
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
1163 # integers. The argument ``check=False`` is required
1164 # because the trace inner-product method for this
1165 # class is a stub and can't actually be checked.
1166 J
= MatrixEuclideanJordanAlgebra(QQ
,
1168 normalize_basis
=False,
1170 a
= J
._charpoly
_coefficients
()
1172 # Unfortunately, changing the basis does change the
1173 # coefficients of the characteristic polynomial, but since
1174 # these are really the coefficients of the "characteristic
1175 # polynomial of" function, everything is still nice and
1176 # unevaluated. It's therefore "obvious" how scaling the
1177 # basis affects the coordinate variables X1, X2, et
1178 # cetera. Scaling the first basis vector up by "n" adds a
1179 # factor of 1/n into every "X1" term, for example. So here
1180 # we simply undo the basis_normalizer scaling that we
1181 # performed earlier.
1183 # The a[0] access here is safe because trivial algebras
1184 # won't have any basis normalizers and therefore won't
1185 # make it to this "else" branch.
1186 XS
= a
[0].parent().gens()
1187 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1188 for i
in range(len(XS
)) }
1189 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1193 def multiplication_table_from_matrix_basis(basis
):
1195 At least three of the five simple Euclidean Jordan algebras have the
1196 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1197 multiplication on the right is matrix multiplication. Given a basis
1198 for the underlying matrix space, this function returns a
1199 multiplication table (obtained by looping through the basis
1200 elements) for an algebra of those matrices.
1202 # In S^2, for example, we nominally have four coordinates even
1203 # though the space is of dimension three only. The vector space V
1204 # is supposed to hold the entire long vector, and the subspace W
1205 # of V will be spanned by the vectors that arise from symmetric
1206 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1210 field
= basis
[0].base_ring()
1211 dimension
= basis
[0].nrows()
1213 V
= VectorSpace(field
, dimension
**2)
1214 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1216 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1219 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1220 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1228 Embed the matrix ``M`` into a space of real matrices.
1230 The matrix ``M`` can have entries in any field at the moment:
1231 the real numbers, complex numbers, or quaternions. And although
1232 they are not a field, we can probably support octonions at some
1233 point, too. This function returns a real matrix that "acts like"
1234 the original with respect to matrix multiplication; i.e.
1236 real_embed(M*N) = real_embed(M)*real_embed(N)
1239 raise NotImplementedError
1243 def real_unembed(M
):
1245 The inverse of :meth:`real_embed`.
1247 raise NotImplementedError
1251 def natural_inner_product(cls
,X
,Y
):
1252 Xu
= cls
.real_unembed(X
)
1253 Yu
= cls
.real_unembed(Y
)
1254 tr
= (Xu
*Yu
).trace()
1257 # It's real already.
1260 # Otherwise, try the thing that works for complex numbers; and
1261 # if that doesn't work, the thing that works for quaternions.
1263 return tr
.vector()[0] # real part, imag part is index 1
1264 except AttributeError:
1265 # A quaternions doesn't have a vector() method, but does
1266 # have coefficient_tuple() method that returns the
1267 # coefficients of 1, i, j, and k -- in that order.
1268 return tr
.coefficient_tuple()[0]
1271 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1275 The identity function, for embedding real matrices into real
1281 def real_unembed(M
):
1283 The identity function, for unembedding real matrices from real
1289 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1291 The rank-n simple EJA consisting of real symmetric n-by-n
1292 matrices, the usual symmetric Jordan product, and the trace inner
1293 product. It has dimension `(n^2 + n)/2` over the reals.
1297 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1301 sage: J = RealSymmetricEJA(2)
1302 sage: e0, e1, e2 = J.gens()
1310 In theory, our "field" can be any subfield of the reals::
1312 sage: RealSymmetricEJA(2, RDF)
1313 Euclidean Jordan algebra of dimension 3 over Real Double Field
1314 sage: RealSymmetricEJA(2, RR)
1315 Euclidean Jordan algebra of dimension 3 over Real Field with
1316 53 bits of precision
1320 The dimension of this algebra is `(n^2 + n) / 2`::
1322 sage: set_random_seed()
1323 sage: n_max = RealSymmetricEJA._max_test_case_size()
1324 sage: n = ZZ.random_element(1, n_max)
1325 sage: J = RealSymmetricEJA(n)
1326 sage: J.dimension() == (n^2 + n)/2
1329 The Jordan multiplication is what we think it is::
1331 sage: set_random_seed()
1332 sage: J = RealSymmetricEJA.random_instance()
1333 sage: x,y = J.random_elements(2)
1334 sage: actual = (x*y).natural_representation()
1335 sage: X = x.natural_representation()
1336 sage: Y = y.natural_representation()
1337 sage: expected = (X*Y + Y*X)/2
1338 sage: actual == expected
1340 sage: J(expected) == x*y
1343 We can change the generator prefix::
1345 sage: RealSymmetricEJA(3, prefix='q').gens()
1346 (q0, q1, q2, q3, q4, q5)
1348 Our natural basis is normalized with respect to the natural inner
1349 product unless we specify otherwise::
1351 sage: set_random_seed()
1352 sage: J = RealSymmetricEJA.random_instance()
1353 sage: all( b.norm() == 1 for b in J.gens() )
1356 Since our natural basis is normalized with respect to the natural
1357 inner product, and since we know that this algebra is an EJA, any
1358 left-multiplication operator's matrix will be symmetric because
1359 natural->EJA basis representation is an isometry and within the EJA
1360 the operator is self-adjoint by the Jordan axiom::
1362 sage: set_random_seed()
1363 sage: x = RealSymmetricEJA.random_instance().random_element()
1364 sage: x.operator().matrix().is_symmetric()
1367 We can construct the (trivial) algebra of rank zero::
1369 sage: RealSymmetricEJA(0)
1370 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1374 def _denormalized_basis(cls
, n
, field
):
1376 Return a basis for the space of real symmetric n-by-n matrices.
1380 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1384 sage: set_random_seed()
1385 sage: n = ZZ.random_element(1,5)
1386 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1387 sage: all( M.is_symmetric() for M in B)
1391 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1395 for j
in range(i
+1):
1396 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1400 Sij
= Eij
+ Eij
.transpose()
1406 def _max_test_case_size():
1407 return 4 # Dimension 10
1410 def __init__(self
, n
, field
=AA
, **kwargs
):
1411 basis
= self
._denormalized
_basis
(n
, field
)
1412 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, **kwargs
)
1413 self
.rank
.set_cache(n
)
1416 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1420 Embed the n-by-n complex matrix ``M`` into the space of real
1421 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1422 bi` to the block matrix ``[[a,b],[-b,a]]``.
1426 sage: from mjo.eja.eja_algebra import \
1427 ....: ComplexMatrixEuclideanJordanAlgebra
1431 sage: F = QuadraticField(-1, 'I')
1432 sage: x1 = F(4 - 2*i)
1433 sage: x2 = F(1 + 2*i)
1436 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1437 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1446 Embedding is a homomorphism (isomorphism, in fact)::
1448 sage: set_random_seed()
1449 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1450 sage: n = ZZ.random_element(n_max)
1451 sage: F = QuadraticField(-1, 'I')
1452 sage: X = random_matrix(F, n)
1453 sage: Y = random_matrix(F, n)
1454 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1455 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1456 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1463 raise ValueError("the matrix 'M' must be square")
1465 # We don't need any adjoined elements...
1466 field
= M
.base_ring().base_ring()
1470 a
= z
.list()[0] # real part, I guess
1471 b
= z
.list()[1] # imag part, I guess
1472 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1474 return matrix
.block(field
, n
, blocks
)
1478 def real_unembed(M
):
1480 The inverse of _embed_complex_matrix().
1484 sage: from mjo.eja.eja_algebra import \
1485 ....: ComplexMatrixEuclideanJordanAlgebra
1489 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1490 ....: [-2, 1, -4, 3],
1491 ....: [ 9, 10, 11, 12],
1492 ....: [-10, 9, -12, 11] ])
1493 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1495 [ 10*I + 9 12*I + 11]
1499 Unembedding is the inverse of embedding::
1501 sage: set_random_seed()
1502 sage: F = QuadraticField(-1, 'I')
1503 sage: M = random_matrix(F, 3)
1504 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1505 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1511 raise ValueError("the matrix 'M' must be square")
1512 if not n
.mod(2).is_zero():
1513 raise ValueError("the matrix 'M' must be a complex embedding")
1515 # If "M" was normalized, its base ring might have roots
1516 # adjoined and they can stick around after unembedding.
1517 field
= M
.base_ring()
1518 R
= PolynomialRing(field
, 'z')
1521 # Sage doesn't know how to embed AA into QQbar, i.e. how
1522 # to adjoin sqrt(-1) to AA.
1525 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1528 # Go top-left to bottom-right (reading order), converting every
1529 # 2-by-2 block we see to a single complex element.
1531 for k
in range(n
/2):
1532 for j
in range(n
/2):
1533 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1534 if submat
[0,0] != submat
[1,1]:
1535 raise ValueError('bad on-diagonal submatrix')
1536 if submat
[0,1] != -submat
[1,0]:
1537 raise ValueError('bad off-diagonal submatrix')
1538 z
= submat
[0,0] + submat
[0,1]*i
1541 return matrix(F
, n
/2, elements
)
1545 def natural_inner_product(cls
,X
,Y
):
1547 Compute a natural inner product in this algebra directly from
1552 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1556 This gives the same answer as the slow, default method implemented
1557 in :class:`MatrixEuclideanJordanAlgebra`::
1559 sage: set_random_seed()
1560 sage: J = ComplexHermitianEJA.random_instance()
1561 sage: x,y = J.random_elements(2)
1562 sage: Xe = x.natural_representation()
1563 sage: Ye = y.natural_representation()
1564 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1565 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1566 sage: expected = (X*Y).trace().real()
1567 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1568 sage: actual == expected
1572 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1575 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1577 The rank-n simple EJA consisting of complex Hermitian n-by-n
1578 matrices over the real numbers, the usual symmetric Jordan product,
1579 and the real-part-of-trace inner product. It has dimension `n^2` over
1584 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1588 In theory, our "field" can be any subfield of the reals::
1590 sage: ComplexHermitianEJA(2, RDF)
1591 Euclidean Jordan algebra of dimension 4 over Real Double Field
1592 sage: ComplexHermitianEJA(2, RR)
1593 Euclidean Jordan algebra of dimension 4 over Real Field with
1594 53 bits of precision
1598 The dimension of this algebra is `n^2`::
1600 sage: set_random_seed()
1601 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1602 sage: n = ZZ.random_element(1, n_max)
1603 sage: J = ComplexHermitianEJA(n)
1604 sage: J.dimension() == n^2
1607 The Jordan multiplication is what we think it is::
1609 sage: set_random_seed()
1610 sage: J = ComplexHermitianEJA.random_instance()
1611 sage: x,y = J.random_elements(2)
1612 sage: actual = (x*y).natural_representation()
1613 sage: X = x.natural_representation()
1614 sage: Y = y.natural_representation()
1615 sage: expected = (X*Y + Y*X)/2
1616 sage: actual == expected
1618 sage: J(expected) == x*y
1621 We can change the generator prefix::
1623 sage: ComplexHermitianEJA(2, prefix='z').gens()
1626 Our natural basis is normalized with respect to the natural inner
1627 product unless we specify otherwise::
1629 sage: set_random_seed()
1630 sage: J = ComplexHermitianEJA.random_instance()
1631 sage: all( b.norm() == 1 for b in J.gens() )
1634 Since our natural basis is normalized with respect to the natural
1635 inner product, and since we know that this algebra is an EJA, any
1636 left-multiplication operator's matrix will be symmetric because
1637 natural->EJA basis representation is an isometry and within the EJA
1638 the operator is self-adjoint by the Jordan axiom::
1640 sage: set_random_seed()
1641 sage: x = ComplexHermitianEJA.random_instance().random_element()
1642 sage: x.operator().matrix().is_symmetric()
1645 We can construct the (trivial) algebra of rank zero::
1647 sage: ComplexHermitianEJA(0)
1648 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1653 def _denormalized_basis(cls
, n
, field
):
1655 Returns a basis for the space of complex Hermitian n-by-n matrices.
1657 Why do we embed these? Basically, because all of numerical linear
1658 algebra assumes that you're working with vectors consisting of `n`
1659 entries from a field and scalars from the same field. There's no way
1660 to tell SageMath that (for example) the vectors contain complex
1661 numbers, while the scalar field is real.
1665 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1669 sage: set_random_seed()
1670 sage: n = ZZ.random_element(1,5)
1671 sage: field = QuadraticField(2, 'sqrt2')
1672 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1673 sage: all( M.is_symmetric() for M in B)
1677 R
= PolynomialRing(field
, 'z')
1679 F
= field
.extension(z
**2 + 1, 'I')
1682 # This is like the symmetric case, but we need to be careful:
1684 # * We want conjugate-symmetry, not just symmetry.
1685 # * The diagonal will (as a result) be real.
1689 for j
in range(i
+1):
1690 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1692 Sij
= cls
.real_embed(Eij
)
1695 # The second one has a minus because it's conjugated.
1696 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1698 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1701 # Since we embedded these, we can drop back to the "field" that we
1702 # started with instead of the complex extension "F".
1703 return ( s
.change_ring(field
) for s
in S
)
1706 def __init__(self
, n
, field
=AA
, **kwargs
):
1707 basis
= self
._denormalized
_basis
(n
,field
)
1708 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1709 self
.rank
.set_cache(n
)
1712 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1716 Embed the n-by-n quaternion matrix ``M`` into the space of real
1717 matrices of size 4n-by-4n by first sending each quaternion entry `z
1718 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1719 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1724 sage: from mjo.eja.eja_algebra import \
1725 ....: QuaternionMatrixEuclideanJordanAlgebra
1729 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1730 sage: i,j,k = Q.gens()
1731 sage: x = 1 + 2*i + 3*j + 4*k
1732 sage: M = matrix(Q, 1, [[x]])
1733 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1739 Embedding is a homomorphism (isomorphism, in fact)::
1741 sage: set_random_seed()
1742 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1743 sage: n = ZZ.random_element(n_max)
1744 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1745 sage: X = random_matrix(Q, n)
1746 sage: Y = random_matrix(Q, n)
1747 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1748 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1749 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1754 quaternions
= M
.base_ring()
1757 raise ValueError("the matrix 'M' must be square")
1759 F
= QuadraticField(-1, 'I')
1764 t
= z
.coefficient_tuple()
1769 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1770 [-c
+ d
*i
, a
- b
*i
]])
1771 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1772 blocks
.append(realM
)
1774 # We should have real entries by now, so use the realest field
1775 # we've got for the return value.
1776 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1781 def real_unembed(M
):
1783 The inverse of _embed_quaternion_matrix().
1787 sage: from mjo.eja.eja_algebra import \
1788 ....: QuaternionMatrixEuclideanJordanAlgebra
1792 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1793 ....: [-2, 1, -4, 3],
1794 ....: [-3, 4, 1, -2],
1795 ....: [-4, -3, 2, 1]])
1796 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1797 [1 + 2*i + 3*j + 4*k]
1801 Unembedding is the inverse of embedding::
1803 sage: set_random_seed()
1804 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1805 sage: M = random_matrix(Q, 3)
1806 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1807 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1813 raise ValueError("the matrix 'M' must be square")
1814 if not n
.mod(4).is_zero():
1815 raise ValueError("the matrix 'M' must be a quaternion embedding")
1817 # Use the base ring of the matrix to ensure that its entries can be
1818 # multiplied by elements of the quaternion algebra.
1819 field
= M
.base_ring()
1820 Q
= QuaternionAlgebra(field
,-1,-1)
1823 # Go top-left to bottom-right (reading order), converting every
1824 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1827 for l
in range(n
/4):
1828 for m
in range(n
/4):
1829 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1830 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1831 if submat
[0,0] != submat
[1,1].conjugate():
1832 raise ValueError('bad on-diagonal submatrix')
1833 if submat
[0,1] != -submat
[1,0].conjugate():
1834 raise ValueError('bad off-diagonal submatrix')
1835 z
= submat
[0,0].real()
1836 z
+= submat
[0,0].imag()*i
1837 z
+= submat
[0,1].real()*j
1838 z
+= submat
[0,1].imag()*k
1841 return matrix(Q
, n
/4, elements
)
1845 def natural_inner_product(cls
,X
,Y
):
1847 Compute a natural inner product in this algebra directly from
1852 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1856 This gives the same answer as the slow, default method implemented
1857 in :class:`MatrixEuclideanJordanAlgebra`::
1859 sage: set_random_seed()
1860 sage: J = QuaternionHermitianEJA.random_instance()
1861 sage: x,y = J.random_elements(2)
1862 sage: Xe = x.natural_representation()
1863 sage: Ye = y.natural_representation()
1864 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1865 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1866 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1867 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1868 sage: actual == expected
1872 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1875 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1877 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1878 matrices, the usual symmetric Jordan product, and the
1879 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1884 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1888 In theory, our "field" can be any subfield of the reals::
1890 sage: QuaternionHermitianEJA(2, RDF)
1891 Euclidean Jordan algebra of dimension 6 over Real Double Field
1892 sage: QuaternionHermitianEJA(2, RR)
1893 Euclidean Jordan algebra of dimension 6 over Real Field with
1894 53 bits of precision
1898 The dimension of this algebra is `2*n^2 - n`::
1900 sage: set_random_seed()
1901 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1902 sage: n = ZZ.random_element(1, n_max)
1903 sage: J = QuaternionHermitianEJA(n)
1904 sage: J.dimension() == 2*(n^2) - n
1907 The Jordan multiplication is what we think it is::
1909 sage: set_random_seed()
1910 sage: J = QuaternionHermitianEJA.random_instance()
1911 sage: x,y = J.random_elements(2)
1912 sage: actual = (x*y).natural_representation()
1913 sage: X = x.natural_representation()
1914 sage: Y = y.natural_representation()
1915 sage: expected = (X*Y + Y*X)/2
1916 sage: actual == expected
1918 sage: J(expected) == x*y
1921 We can change the generator prefix::
1923 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1924 (a0, a1, a2, a3, a4, a5)
1926 Our natural basis is normalized with respect to the natural inner
1927 product unless we specify otherwise::
1929 sage: set_random_seed()
1930 sage: J = QuaternionHermitianEJA.random_instance()
1931 sage: all( b.norm() == 1 for b in J.gens() )
1934 Since our natural basis is normalized with respect to the natural
1935 inner product, and since we know that this algebra is an EJA, any
1936 left-multiplication operator's matrix will be symmetric because
1937 natural->EJA basis representation is an isometry and within the EJA
1938 the operator is self-adjoint by the Jordan axiom::
1940 sage: set_random_seed()
1941 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1942 sage: x.operator().matrix().is_symmetric()
1945 We can construct the (trivial) algebra of rank zero::
1947 sage: QuaternionHermitianEJA(0)
1948 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1952 def _denormalized_basis(cls
, n
, field
):
1954 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1956 Why do we embed these? Basically, because all of numerical
1957 linear algebra assumes that you're working with vectors consisting
1958 of `n` entries from a field and scalars from the same field. There's
1959 no way to tell SageMath that (for example) the vectors contain
1960 complex numbers, while the scalar field is real.
1964 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1968 sage: set_random_seed()
1969 sage: n = ZZ.random_element(1,5)
1970 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1971 sage: all( M.is_symmetric() for M in B )
1975 Q
= QuaternionAlgebra(QQ
,-1,-1)
1978 # This is like the symmetric case, but we need to be careful:
1980 # * We want conjugate-symmetry, not just symmetry.
1981 # * The diagonal will (as a result) be real.
1985 for j
in range(i
+1):
1986 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1988 Sij
= cls
.real_embed(Eij
)
1991 # The second, third, and fourth ones have a minus
1992 # because they're conjugated.
1993 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1995 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1997 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1999 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2002 # Since we embedded these, we can drop back to the "field" that we
2003 # started with instead of the quaternion algebra "Q".
2004 return ( s
.change_ring(field
) for s
in S
)
2007 def __init__(self
, n
, field
=AA
, **kwargs
):
2008 basis
= self
._denormalized
_basis
(n
,field
)
2009 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
2010 self
.rank
.set_cache(n
)
2013 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2015 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2016 with the half-trace inner product and jordan product ``x*y =
2017 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
2018 symmetric positive-definite "bilinear form" matrix. It has
2019 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
2020 when ``B`` is the identity matrix of order ``n-1``.
2024 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2025 ....: JordanSpinEJA)
2029 When no bilinear form is specified, the identity matrix is used,
2030 and the resulting algebra is the Jordan spin algebra::
2032 sage: J0 = BilinearFormEJA(3)
2033 sage: J1 = JordanSpinEJA(3)
2034 sage: J0.multiplication_table() == J0.multiplication_table()
2039 We can create a zero-dimensional algebra::
2041 sage: J = BilinearFormEJA(0)
2045 We can check the multiplication condition given in the Jordan, von
2046 Neumann, and Wigner paper (and also discussed on my "On the
2047 symmetry..." paper). Note that this relies heavily on the standard
2048 choice of basis, as does anything utilizing the bilinear form matrix::
2050 sage: set_random_seed()
2051 sage: n = ZZ.random_element(5)
2052 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2053 sage: B = M.transpose()*M
2054 sage: J = BilinearFormEJA(n, B=B)
2055 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2056 sage: V = J.vector_space()
2057 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2058 ....: for ei in eis ]
2059 sage: actual = [ sis[i]*sis[j]
2060 ....: for i in range(n-1)
2061 ....: for j in range(n-1) ]
2062 sage: expected = [ J.one() if i == j else J.zero()
2063 ....: for i in range(n-1)
2064 ....: for j in range(n-1) ]
2065 sage: actual == expected
2068 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
2070 self
._B
= matrix
.identity(field
, max(0,n
-1))
2074 V
= VectorSpace(field
, n
)
2075 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2084 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
2085 zbar
= y0
*xbar
+ x0
*ybar
2086 z
= V([z0
] + zbar
.list())
2087 mult_table
[i
][j
] = z
2089 # The rank of this algebra is two, unless we're in a
2090 # one-dimensional ambient space (because the rank is bounded
2091 # by the ambient dimension).
2092 fdeja
= super(BilinearFormEJA
, self
)
2093 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2094 self
.rank
.set_cache(min(n
,2))
2096 def inner_product(self
, x
, y
):
2098 Half of the trace inner product.
2100 This is defined so that the special case of the Jordan spin
2101 algebra gets the usual inner product.
2105 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2109 Ensure that this is one-half of the trace inner-product when
2110 the algebra isn't just the reals (when ``n`` isn't one). This
2111 is in Faraut and Koranyi, and also my "On the symmetry..."
2114 sage: set_random_seed()
2115 sage: n = ZZ.random_element(2,5)
2116 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2117 sage: B = M.transpose()*M
2118 sage: J = BilinearFormEJA(n, B=B)
2119 sage: x = J.random_element()
2120 sage: y = J.random_element()
2121 sage: x.inner_product(y) == (x*y).trace()/2
2125 xvec
= x
.to_vector()
2127 yvec
= y
.to_vector()
2129 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
2132 class JordanSpinEJA(BilinearFormEJA
):
2134 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2135 with the usual inner product and jordan product ``x*y =
2136 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2141 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2145 This multiplication table can be verified by hand::
2147 sage: J = JordanSpinEJA(4)
2148 sage: e0,e1,e2,e3 = J.gens()
2164 We can change the generator prefix::
2166 sage: JordanSpinEJA(2, prefix='B').gens()
2171 Ensure that we have the usual inner product on `R^n`::
2173 sage: set_random_seed()
2174 sage: J = JordanSpinEJA.random_instance()
2175 sage: x,y = J.random_elements(2)
2176 sage: X = x.natural_representation()
2177 sage: Y = y.natural_representation()
2178 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2182 def __init__(self
, n
, field
=AA
, **kwargs
):
2183 # This is a special case of the BilinearFormEJA with the identity
2184 # matrix as its bilinear form.
2185 return super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2188 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2190 The trivial Euclidean Jordan algebra consisting of only a zero element.
2194 sage: from mjo.eja.eja_algebra import TrivialEJA
2198 sage: J = TrivialEJA()
2205 sage: 7*J.one()*12*J.one()
2207 sage: J.one().inner_product(J.one())
2209 sage: J.one().norm()
2211 sage: J.one().subalgebra_generated_by()
2212 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2217 def __init__(self
, field
=AA
, **kwargs
):
2219 fdeja
= super(TrivialEJA
, self
)
2220 # The rank is zero using my definition, namely the dimension of the
2221 # largest subalgebra generated by any element.
2222 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2223 self
.rank
.set_cache(0)
2226 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2228 The external (orthogonal) direct sum of two other Euclidean Jordan
2229 algebras. Essentially the Cartesian product of its two factors.
2230 Every Euclidean Jordan algebra decomposes into an orthogonal
2231 direct sum of simple Euclidean Jordan algebras, so no generality
2232 is lost by providing only this construction.
2236 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2237 ....: RealSymmetricEJA,
2242 sage: J1 = HadamardEJA(2)
2243 sage: J2 = RealSymmetricEJA(3)
2244 sage: J = DirectSumEJA(J1,J2)
2251 def __init__(self
, J1
, J2
, field
=AA
, **kwargs
):
2255 V
= VectorSpace(field
, n
)
2256 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2260 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2261 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2265 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2266 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2268 fdeja
= super(DirectSumEJA
, self
)
2269 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2270 self
.rank
.set_cache(J1
.rank() + J2
.rank())