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.
299 for i
in range(self
.dimension()):
300 for j
in range(self
.dimension()):
301 for k
in range(self
.dimension()):
305 if (x
*y
).inner_product(z
) != x
.inner_product(y
*z
):
311 def characteristic_polynomial_of(self
):
313 Return the algebra's "characteristic polynomial of" function,
314 which is itself a multivariate polynomial that, when evaluated
315 at the coordinates of some algebra element, returns that
316 element's characteristic polynomial.
318 The resulting polynomial has `n+1` variables, where `n` is the
319 dimension of this algebra. The first `n` variables correspond to
320 the coordinates of an algebra element: when evaluated at the
321 coordinates of an algebra element with respect to a certain
322 basis, the result is a univariate polynomial (in the one
323 remaining variable ``t``), namely the characteristic polynomial
328 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
332 The characteristic polynomial in the spin algebra is given in
333 Alizadeh, Example 11.11::
335 sage: J = JordanSpinEJA(3)
336 sage: p = J.characteristic_polynomial_of(); p
337 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
338 sage: xvec = J.one().to_vector()
342 By definition, the characteristic polynomial is a monic
343 degree-zero polynomial in a rank-zero algebra. Note that
344 Cayley-Hamilton is indeed satisfied since the polynomial
345 ``1`` evaluates to the identity element of the algebra on
348 sage: J = TrivialEJA()
349 sage: J.characteristic_polynomial_of()
356 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
357 a
= self
._charpoly
_coefficients
()
359 # We go to a bit of trouble here to reorder the
360 # indeterminates, so that it's easier to evaluate the
361 # characteristic polynomial at x's coordinates and get back
362 # something in terms of t, which is what we want.
363 S
= PolynomialRing(self
.base_ring(),'t')
367 S
= PolynomialRing(S
, R
.variable_names())
370 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
373 def inner_product(self
, x
, y
):
375 The inner product associated with this Euclidean Jordan algebra.
377 Defaults to the trace inner product, but can be overridden by
378 subclasses if they are sure that the necessary properties are
383 sage: from mjo.eja.eja_algebra import random_eja
387 Our inner product is "associative," which means the following for
388 a symmetric bilinear form::
390 sage: set_random_seed()
391 sage: J = random_eja()
392 sage: x,y,z = J.random_elements(3)
393 sage: (x*y).inner_product(z) == y.inner_product(x*z)
397 X
= x
.natural_representation()
398 Y
= y
.natural_representation()
399 return self
.natural_inner_product(X
,Y
)
402 def is_trivial(self
):
404 Return whether or not this algebra is trivial.
406 A trivial algebra contains only the zero element.
410 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
415 sage: J = ComplexHermitianEJA(3)
421 sage: J = TrivialEJA()
426 return self
.dimension() == 0
429 def multiplication_table(self
):
431 Return a visual representation of this algebra's multiplication
432 table (on basis elements).
436 sage: from mjo.eja.eja_algebra import JordanSpinEJA
440 sage: J = JordanSpinEJA(4)
441 sage: J.multiplication_table()
442 +----++----+----+----+----+
443 | * || e0 | e1 | e2 | e3 |
444 +====++====+====+====+====+
445 | e0 || e0 | e1 | e2 | e3 |
446 +----++----+----+----+----+
447 | e1 || e1 | e0 | 0 | 0 |
448 +----++----+----+----+----+
449 | e2 || e2 | 0 | e0 | 0 |
450 +----++----+----+----+----+
451 | e3 || e3 | 0 | 0 | e0 |
452 +----++----+----+----+----+
455 M
= list(self
._multiplication
_table
) # copy
456 for i
in range(len(M
)):
457 # M had better be "square"
458 M
[i
] = [self
.monomial(i
)] + M
[i
]
459 M
= [["*"] + list(self
.gens())] + M
460 return table(M
, header_row
=True, header_column
=True, frame
=True)
463 def natural_basis(self
):
465 Return a more-natural representation of this algebra's basis.
467 Every finite-dimensional Euclidean Jordan Algebra is a direct
468 sum of five simple algebras, four of which comprise Hermitian
469 matrices. This method returns the original "natural" basis
470 for our underlying vector space. (Typically, the natural basis
471 is used to construct the multiplication table in the first place.)
473 Note that this will always return a matrix. The standard basis
474 in `R^n` will be returned as `n`-by-`1` column matrices.
478 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
479 ....: RealSymmetricEJA)
483 sage: J = RealSymmetricEJA(2)
485 Finite family {0: e0, 1: e1, 2: e2}
486 sage: J.natural_basis()
488 [1 0] [ 0 0.7071067811865475?] [0 0]
489 [0 0], [0.7071067811865475? 0], [0 1]
494 sage: J = JordanSpinEJA(2)
496 Finite family {0: e0, 1: e1}
497 sage: J.natural_basis()
504 if self
._natural
_basis
is None:
505 M
= self
.natural_basis_space()
506 return tuple( M(b
.to_vector()) for b
in self
.basis() )
508 return self
._natural
_basis
511 def natural_basis_space(self
):
513 Return the matrix space in which this algebra's natural basis
516 Generally this will be an `n`-by-`1` column-vector space,
517 except when the algebra is trivial. There it's `n`-by-`n`
518 (where `n` is zero), to ensure that two elements of the
519 natural basis space (empty matrices) can be multiplied.
521 if self
.is_trivial():
522 return MatrixSpace(self
.base_ring(), 0)
523 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
524 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
526 return self
._natural
_basis
[0].matrix_space()
530 def natural_inner_product(X
,Y
):
532 Compute the inner product of two naturally-represented elements.
534 For example in the real symmetric matrix EJA, this will compute
535 the trace inner-product of two n-by-n symmetric matrices. The
536 default should work for the real cartesian product EJA, the
537 Jordan spin EJA, and the real symmetric matrices. The others
538 will have to be overridden.
540 return (X
.conjugate_transpose()*Y
).trace()
546 Return the unit element of this algebra.
550 sage: from mjo.eja.eja_algebra import (HadamardEJA,
555 sage: J = HadamardEJA(5)
557 e0 + e1 + e2 + e3 + e4
561 The identity element acts like the identity::
563 sage: set_random_seed()
564 sage: J = random_eja()
565 sage: x = J.random_element()
566 sage: J.one()*x == x and x*J.one() == x
569 The matrix of the unit element's operator is the identity::
571 sage: set_random_seed()
572 sage: J = random_eja()
573 sage: actual = J.one().operator().matrix()
574 sage: expected = matrix.identity(J.base_ring(), J.dimension())
575 sage: actual == expected
579 # We can brute-force compute the matrices of the operators
580 # that correspond to the basis elements of this algebra.
581 # If some linear combination of those basis elements is the
582 # algebra identity, then the same linear combination of
583 # their matrices has to be the identity matrix.
585 # Of course, matrices aren't vectors in sage, so we have to
586 # appeal to the "long vectors" isometry.
587 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
589 # Now we use basis linear algebra to find the coefficients,
590 # of the matrices-as-vectors-linear-combination, which should
591 # work for the original algebra basis too.
592 A
= matrix
.column(self
.base_ring(), oper_vecs
)
594 # We used the isometry on the left-hand side already, but we
595 # still need to do it for the right-hand side. Recall that we
596 # wanted something that summed to the identity matrix.
597 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
599 # Now if there's an identity element in the algebra, this should work.
600 coeffs
= A
.solve_right(b
)
601 return self
.linear_combination(zip(self
.gens(), coeffs
))
604 def peirce_decomposition(self
, c
):
606 The Peirce decomposition of this algebra relative to the
609 In the future, this can be extended to a complete system of
610 orthogonal idempotents.
614 - ``c`` -- an idempotent of this algebra.
618 A triple (J0, J5, J1) containing two subalgebras and one subspace
621 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
622 corresponding to the eigenvalue zero.
624 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
625 corresponding to the eigenvalue one-half.
627 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
628 corresponding to the eigenvalue one.
630 These are the only possible eigenspaces for that operator, and this
631 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
632 orthogonal, and are subalgebras of this algebra with the appropriate
637 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
641 The canonical example comes from the symmetric matrices, which
642 decompose into diagonal and off-diagonal parts::
644 sage: J = RealSymmetricEJA(3)
645 sage: C = matrix(QQ, [ [1,0,0],
649 sage: J0,J5,J1 = J.peirce_decomposition(c)
651 Euclidean Jordan algebra of dimension 1...
653 Vector space of degree 6 and dimension 2...
655 Euclidean Jordan algebra of dimension 3...
656 sage: J0.one().natural_representation()
660 sage: orig_df = AA.options.display_format
661 sage: AA.options.display_format = 'radical'
662 sage: J.from_vector(J5.basis()[0]).natural_representation()
666 sage: J.from_vector(J5.basis()[1]).natural_representation()
670 sage: AA.options.display_format = orig_df
671 sage: J1.one().natural_representation()
678 Every algebra decomposes trivially with respect to its identity
681 sage: set_random_seed()
682 sage: J = random_eja()
683 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
684 sage: J0.dimension() == 0 and J5.dimension() == 0
686 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
689 The decomposition is into eigenspaces, and its components are
690 therefore necessarily orthogonal. Moreover, the identity
691 elements in the two subalgebras are the projections onto their
692 respective subspaces of the superalgebra's identity element::
694 sage: set_random_seed()
695 sage: J = random_eja()
696 sage: x = J.random_element()
697 sage: if not J.is_trivial():
698 ....: while x.is_nilpotent():
699 ....: x = J.random_element()
700 sage: c = x.subalgebra_idempotent()
701 sage: J0,J5,J1 = J.peirce_decomposition(c)
703 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
704 ....: w = w.superalgebra_element()
705 ....: y = J.from_vector(y)
706 ....: z = z.superalgebra_element()
707 ....: ipsum += w.inner_product(y).abs()
708 ....: ipsum += w.inner_product(z).abs()
709 ....: ipsum += y.inner_product(z).abs()
712 sage: J1(c) == J1.one()
714 sage: J0(J.one() - c) == J0.one()
718 if not c
.is_idempotent():
719 raise ValueError("element is not idempotent: %s" % c
)
721 # Default these to what they should be if they turn out to be
722 # trivial, because eigenspaces_left() won't return eigenvalues
723 # corresponding to trivial spaces (e.g. it returns only the
724 # eigenspace corresponding to lambda=1 if you take the
725 # decomposition relative to the identity element).
726 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
727 J0
= trivial
# eigenvalue zero
728 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
729 J1
= trivial
# eigenvalue one
731 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
732 if eigval
== ~
(self
.base_ring()(2)):
735 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
736 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
742 raise ValueError("unexpected eigenvalue: %s" % eigval
)
747 def random_element(self
, thorough
=False):
749 Return a random element of this algebra.
751 Our algebra superclass method only returns a linear
752 combination of at most two basis elements. We instead
753 want the vector space "random element" method that
754 returns a more diverse selection.
758 - ``thorough`` -- (boolean; default False) whether or not we
759 should generate irrational coefficients for the random
760 element when our base ring is irrational; this slows the
761 algebra operations to a crawl, but any truly random method
765 # For a general base ring... maybe we can trust this to do the
766 # right thing? Unlikely, but.
767 V
= self
.vector_space()
768 v
= V
.random_element()
770 if self
.base_ring() is AA
:
771 # The "random element" method of the algebraic reals is
772 # stupid at the moment, and only returns integers between
773 # -2 and 2, inclusive:
775 # https://trac.sagemath.org/ticket/30875
777 # Instead, we implement our own "random vector" method,
778 # and then coerce that into the algebra. We use the vector
779 # space degree here instead of the dimension because a
780 # subalgebra could (for example) be spanned by only two
781 # vectors, each with five coordinates. We need to
782 # generate all five coordinates.
784 v
*= QQbar
.random_element().real()
786 v
*= QQ
.random_element()
788 return self
.from_vector(V
.coordinate_vector(v
))
790 def random_elements(self
, count
, thorough
=False):
792 Return ``count`` random elements as a tuple.
796 - ``thorough`` -- (boolean; default False) whether or not we
797 should generate irrational coefficients for the random
798 elements when our base ring is irrational; this slows the
799 algebra operations to a crawl, but any truly random method
804 sage: from mjo.eja.eja_algebra import JordanSpinEJA
808 sage: J = JordanSpinEJA(3)
809 sage: x,y,z = J.random_elements(3)
810 sage: all( [ x in J, y in J, z in J ])
812 sage: len( J.random_elements(10) ) == 10
816 return tuple( self
.random_element(thorough
)
817 for idx
in range(count
) )
820 def random_instance(cls
, field
=AA
, **kwargs
):
822 Return a random instance of this type of algebra.
824 Beware, this will crash for "most instances" because the
825 constructor below looks wrong.
827 if cls
is TrivialEJA
:
828 # The TrivialEJA class doesn't take an "n" argument because
832 n
= ZZ
.random_element(cls
._max
_test
_case
_size
() + 1)
833 return cls(n
, field
, **kwargs
)
836 def _charpoly_coefficients(self
):
838 The `r` polynomial coefficients of the "characteristic polynomial
842 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
843 R
= PolynomialRing(self
.base_ring(), var_names
)
845 F
= R
.fraction_field()
848 # From a result in my book, these are the entries of the
849 # basis representation of L_x.
850 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
853 L_x
= matrix(F
, n
, n
, L_x_i_j
)
856 if self
.rank
.is_in_cache():
858 # There's no need to pad the system with redundant
859 # columns if we *know* they'll be redundant.
862 # Compute an extra power in case the rank is equal to
863 # the dimension (otherwise, we would stop at x^(r-1)).
864 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
865 for k
in range(n
+1) ]
866 A
= matrix
.column(F
, x_powers
[:n
])
867 AE
= A
.extended_echelon_form()
874 # The theory says that only the first "r" coefficients are
875 # nonzero, and they actually live in the original polynomial
876 # ring and not the fraction field. We negate them because
877 # in the actual characteristic polynomial, they get moved
878 # to the other side where x^r lives.
879 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
884 Return the rank of this EJA.
886 This is a cached method because we know the rank a priori for
887 all of the algebras we can construct. Thus we can avoid the
888 expensive ``_charpoly_coefficients()`` call unless we truly
889 need to compute the whole characteristic polynomial.
893 sage: from mjo.eja.eja_algebra import (HadamardEJA,
895 ....: RealSymmetricEJA,
896 ....: ComplexHermitianEJA,
897 ....: QuaternionHermitianEJA,
902 The rank of the Jordan spin algebra is always two::
904 sage: JordanSpinEJA(2).rank()
906 sage: JordanSpinEJA(3).rank()
908 sage: JordanSpinEJA(4).rank()
911 The rank of the `n`-by-`n` Hermitian real, complex, or
912 quaternion matrices is `n`::
914 sage: RealSymmetricEJA(4).rank()
916 sage: ComplexHermitianEJA(3).rank()
918 sage: QuaternionHermitianEJA(2).rank()
923 Ensure that every EJA that we know how to construct has a
924 positive integer rank, unless the algebra is trivial in
925 which case its rank will be zero::
927 sage: set_random_seed()
928 sage: J = random_eja()
932 sage: r > 0 or (r == 0 and J.is_trivial())
935 Ensure that computing the rank actually works, since the ranks
936 of all simple algebras are known and will be cached by default::
938 sage: J = HadamardEJA(4)
939 sage: J.rank.clear_cache()
945 sage: J = JordanSpinEJA(4)
946 sage: J.rank.clear_cache()
952 sage: J = RealSymmetricEJA(3)
953 sage: J.rank.clear_cache()
959 sage: J = ComplexHermitianEJA(2)
960 sage: J.rank.clear_cache()
966 sage: J = QuaternionHermitianEJA(2)
967 sage: J.rank.clear_cache()
971 return len(self
._charpoly
_coefficients
())
974 def vector_space(self
):
976 Return the vector space that underlies this algebra.
980 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
984 sage: J = RealSymmetricEJA(2)
985 sage: J.vector_space()
986 Vector space of dimension 3 over...
989 return self
.zero().to_vector().parent().ambient_vector_space()
992 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
995 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
):
997 Return the Euclidean Jordan Algebra corresponding to the set
998 `R^n` under the Hadamard product.
1000 Note: this is nothing more than the Cartesian product of ``n``
1001 copies of the spin algebra. Once Cartesian product algebras
1002 are implemented, this can go.
1006 sage: from mjo.eja.eja_algebra import HadamardEJA
1010 This multiplication table can be verified by hand::
1012 sage: J = HadamardEJA(3)
1013 sage: e0,e1,e2 = J.gens()
1029 We can change the generator prefix::
1031 sage: HadamardEJA(3, prefix='r').gens()
1035 def __init__(self
, n
, field
=AA
, **kwargs
):
1036 V
= VectorSpace(field
, n
)
1037 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
1040 fdeja
= super(HadamardEJA
, self
)
1041 fdeja
.__init
__(field
, mult_table
, **kwargs
)
1042 self
.rank
.set_cache(n
)
1044 def inner_product(self
, x
, y
):
1046 Faster to reimplement than to use natural representations.
1050 sage: from mjo.eja.eja_algebra import HadamardEJA
1054 Ensure that this is the usual inner product for the algebras
1057 sage: set_random_seed()
1058 sage: J = HadamardEJA.random_instance()
1059 sage: x,y = J.random_elements(2)
1060 sage: X = x.natural_representation()
1061 sage: Y = y.natural_representation()
1062 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1066 return x
.to_vector().inner_product(y
.to_vector())
1069 def random_eja(field
=AA
):
1071 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1075 sage: from mjo.eja.eja_algebra import random_eja
1080 Euclidean Jordan algebra of dimension...
1083 classname
= choice([TrivialEJA
,
1087 ComplexHermitianEJA
,
1088 QuaternionHermitianEJA
])
1089 return classname
.random_instance(field
=field
)
1094 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1096 def _max_test_case_size():
1097 # Play it safe, since this will be squared and the underlying
1098 # field can have dimension 4 (quaternions) too.
1101 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1103 Compared to the superclass constructor, we take a basis instead of
1104 a multiplication table because the latter can be computed in terms
1105 of the former when the product is known (like it is here).
1107 # Used in this class's fast _charpoly_coefficients() override.
1108 self
._basis
_normalizers
= None
1110 # We're going to loop through this a few times, so now's a good
1111 # time to ensure that it isn't a generator expression.
1112 basis
= tuple(basis
)
1114 if len(basis
) > 1 and normalize_basis
:
1115 # We'll need sqrt(2) to normalize the basis, and this
1116 # winds up in the multiplication table, so the whole
1117 # algebra needs to be over the field extension.
1118 R
= PolynomialRing(field
, 'z')
1121 if p
.is_irreducible():
1122 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1123 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1124 self
._basis
_normalizers
= tuple(
1125 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1126 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1128 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1130 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
1131 fdeja
.__init
__(field
, Qs
, natural_basis
=basis
, **kwargs
)
1136 def _charpoly_coefficients(self
):
1138 Override the parent method with something that tries to compute
1139 over a faster (non-extension) field.
1141 if self
._basis
_normalizers
is None:
1142 # We didn't normalize, so assume that the basis we started
1143 # with had entries in a nice field.
1144 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1146 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1147 self
._basis
_normalizers
) )
1149 # Do this over the rationals and convert back at the end.
1150 # Only works because we know the entries of the basis are
1152 J
= MatrixEuclideanJordanAlgebra(QQ
,
1154 normalize_basis
=False)
1155 a
= J
._charpoly
_coefficients
()
1157 # Unfortunately, changing the basis does change the
1158 # coefficients of the characteristic polynomial, but since
1159 # these are really the coefficients of the "characteristic
1160 # polynomial of" function, everything is still nice and
1161 # unevaluated. It's therefore "obvious" how scaling the
1162 # basis affects the coordinate variables X1, X2, et
1163 # cetera. Scaling the first basis vector up by "n" adds a
1164 # factor of 1/n into every "X1" term, for example. So here
1165 # we simply undo the basis_normalizer scaling that we
1166 # performed earlier.
1168 # The a[0] access here is safe because trivial algebras
1169 # won't have any basis normalizers and therefore won't
1170 # make it to this "else" branch.
1171 XS
= a
[0].parent().gens()
1172 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1173 for i
in range(len(XS
)) }
1174 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1178 def multiplication_table_from_matrix_basis(basis
):
1180 At least three of the five simple Euclidean Jordan algebras have the
1181 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1182 multiplication on the right is matrix multiplication. Given a basis
1183 for the underlying matrix space, this function returns a
1184 multiplication table (obtained by looping through the basis
1185 elements) for an algebra of those matrices.
1187 # In S^2, for example, we nominally have four coordinates even
1188 # though the space is of dimension three only. The vector space V
1189 # is supposed to hold the entire long vector, and the subspace W
1190 # of V will be spanned by the vectors that arise from symmetric
1191 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1195 field
= basis
[0].base_ring()
1196 dimension
= basis
[0].nrows()
1198 V
= VectorSpace(field
, dimension
**2)
1199 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1201 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1204 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1205 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1213 Embed the matrix ``M`` into a space of real matrices.
1215 The matrix ``M`` can have entries in any field at the moment:
1216 the real numbers, complex numbers, or quaternions. And although
1217 they are not a field, we can probably support octonions at some
1218 point, too. This function returns a real matrix that "acts like"
1219 the original with respect to matrix multiplication; i.e.
1221 real_embed(M*N) = real_embed(M)*real_embed(N)
1224 raise NotImplementedError
1228 def real_unembed(M
):
1230 The inverse of :meth:`real_embed`.
1232 raise NotImplementedError
1236 def natural_inner_product(cls
,X
,Y
):
1237 Xu
= cls
.real_unembed(X
)
1238 Yu
= cls
.real_unembed(Y
)
1239 tr
= (Xu
*Yu
).trace()
1242 # It's real already.
1245 # Otherwise, try the thing that works for complex numbers; and
1246 # if that doesn't work, the thing that works for quaternions.
1248 return tr
.vector()[0] # real part, imag part is index 1
1249 except AttributeError:
1250 # A quaternions doesn't have a vector() method, but does
1251 # have coefficient_tuple() method that returns the
1252 # coefficients of 1, i, j, and k -- in that order.
1253 return tr
.coefficient_tuple()[0]
1256 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1260 The identity function, for embedding real matrices into real
1266 def real_unembed(M
):
1268 The identity function, for unembedding real matrices from real
1274 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1276 The rank-n simple EJA consisting of real symmetric n-by-n
1277 matrices, the usual symmetric Jordan product, and the trace inner
1278 product. It has dimension `(n^2 + n)/2` over the reals.
1282 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1286 sage: J = RealSymmetricEJA(2)
1287 sage: e0, e1, e2 = J.gens()
1295 In theory, our "field" can be any subfield of the reals::
1297 sage: RealSymmetricEJA(2, RDF)
1298 Euclidean Jordan algebra of dimension 3 over Real Double Field
1299 sage: RealSymmetricEJA(2, RR)
1300 Euclidean Jordan algebra of dimension 3 over Real Field with
1301 53 bits of precision
1305 The dimension of this algebra is `(n^2 + n) / 2`::
1307 sage: set_random_seed()
1308 sage: n_max = RealSymmetricEJA._max_test_case_size()
1309 sage: n = ZZ.random_element(1, n_max)
1310 sage: J = RealSymmetricEJA(n)
1311 sage: J.dimension() == (n^2 + n)/2
1314 The Jordan multiplication is what we think it is::
1316 sage: set_random_seed()
1317 sage: J = RealSymmetricEJA.random_instance()
1318 sage: x,y = J.random_elements(2)
1319 sage: actual = (x*y).natural_representation()
1320 sage: X = x.natural_representation()
1321 sage: Y = y.natural_representation()
1322 sage: expected = (X*Y + Y*X)/2
1323 sage: actual == expected
1325 sage: J(expected) == x*y
1328 We can change the generator prefix::
1330 sage: RealSymmetricEJA(3, prefix='q').gens()
1331 (q0, q1, q2, q3, q4, q5)
1333 Our natural basis is normalized with respect to the natural inner
1334 product unless we specify otherwise::
1336 sage: set_random_seed()
1337 sage: J = RealSymmetricEJA.random_instance()
1338 sage: all( b.norm() == 1 for b in J.gens() )
1341 Since our natural basis is normalized with respect to the natural
1342 inner product, and since we know that this algebra is an EJA, any
1343 left-multiplication operator's matrix will be symmetric because
1344 natural->EJA basis representation is an isometry and within the EJA
1345 the operator is self-adjoint by the Jordan axiom::
1347 sage: set_random_seed()
1348 sage: x = RealSymmetricEJA.random_instance().random_element()
1349 sage: x.operator().matrix().is_symmetric()
1352 We can construct the (trivial) algebra of rank zero::
1354 sage: RealSymmetricEJA(0)
1355 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1359 def _denormalized_basis(cls
, n
, field
):
1361 Return a basis for the space of real symmetric n-by-n matrices.
1365 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1369 sage: set_random_seed()
1370 sage: n = ZZ.random_element(1,5)
1371 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1372 sage: all( M.is_symmetric() for M in B)
1376 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1380 for j
in range(i
+1):
1381 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1385 Sij
= Eij
+ Eij
.transpose()
1391 def _max_test_case_size():
1392 return 4 # Dimension 10
1395 def __init__(self
, n
, field
=AA
, **kwargs
):
1396 basis
= self
._denormalized
_basis
(n
, field
)
1397 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, **kwargs
)
1398 self
.rank
.set_cache(n
)
1401 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1405 Embed the n-by-n complex matrix ``M`` into the space of real
1406 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1407 bi` to the block matrix ``[[a,b],[-b,a]]``.
1411 sage: from mjo.eja.eja_algebra import \
1412 ....: ComplexMatrixEuclideanJordanAlgebra
1416 sage: F = QuadraticField(-1, 'I')
1417 sage: x1 = F(4 - 2*i)
1418 sage: x2 = F(1 + 2*i)
1421 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1422 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1431 Embedding is a homomorphism (isomorphism, in fact)::
1433 sage: set_random_seed()
1434 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1435 sage: n = ZZ.random_element(n_max)
1436 sage: F = QuadraticField(-1, 'I')
1437 sage: X = random_matrix(F, n)
1438 sage: Y = random_matrix(F, n)
1439 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1440 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1441 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1448 raise ValueError("the matrix 'M' must be square")
1450 # We don't need any adjoined elements...
1451 field
= M
.base_ring().base_ring()
1455 a
= z
.list()[0] # real part, I guess
1456 b
= z
.list()[1] # imag part, I guess
1457 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1459 return matrix
.block(field
, n
, blocks
)
1463 def real_unembed(M
):
1465 The inverse of _embed_complex_matrix().
1469 sage: from mjo.eja.eja_algebra import \
1470 ....: ComplexMatrixEuclideanJordanAlgebra
1474 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1475 ....: [-2, 1, -4, 3],
1476 ....: [ 9, 10, 11, 12],
1477 ....: [-10, 9, -12, 11] ])
1478 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1480 [ 10*I + 9 12*I + 11]
1484 Unembedding is the inverse of embedding::
1486 sage: set_random_seed()
1487 sage: F = QuadraticField(-1, 'I')
1488 sage: M = random_matrix(F, 3)
1489 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1490 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1496 raise ValueError("the matrix 'M' must be square")
1497 if not n
.mod(2).is_zero():
1498 raise ValueError("the matrix 'M' must be a complex embedding")
1500 # If "M" was normalized, its base ring might have roots
1501 # adjoined and they can stick around after unembedding.
1502 field
= M
.base_ring()
1503 R
= PolynomialRing(field
, 'z')
1506 # Sage doesn't know how to embed AA into QQbar, i.e. how
1507 # to adjoin sqrt(-1) to AA.
1510 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1513 # Go top-left to bottom-right (reading order), converting every
1514 # 2-by-2 block we see to a single complex element.
1516 for k
in range(n
/2):
1517 for j
in range(n
/2):
1518 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1519 if submat
[0,0] != submat
[1,1]:
1520 raise ValueError('bad on-diagonal submatrix')
1521 if submat
[0,1] != -submat
[1,0]:
1522 raise ValueError('bad off-diagonal submatrix')
1523 z
= submat
[0,0] + submat
[0,1]*i
1526 return matrix(F
, n
/2, elements
)
1530 def natural_inner_product(cls
,X
,Y
):
1532 Compute a natural inner product in this algebra directly from
1537 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1541 This gives the same answer as the slow, default method implemented
1542 in :class:`MatrixEuclideanJordanAlgebra`::
1544 sage: set_random_seed()
1545 sage: J = ComplexHermitianEJA.random_instance()
1546 sage: x,y = J.random_elements(2)
1547 sage: Xe = x.natural_representation()
1548 sage: Ye = y.natural_representation()
1549 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1550 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1551 sage: expected = (X*Y).trace().real()
1552 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1553 sage: actual == expected
1557 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1560 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1562 The rank-n simple EJA consisting of complex Hermitian n-by-n
1563 matrices over the real numbers, the usual symmetric Jordan product,
1564 and the real-part-of-trace inner product. It has dimension `n^2` over
1569 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1573 In theory, our "field" can be any subfield of the reals::
1575 sage: ComplexHermitianEJA(2, RDF)
1576 Euclidean Jordan algebra of dimension 4 over Real Double Field
1577 sage: ComplexHermitianEJA(2, RR)
1578 Euclidean Jordan algebra of dimension 4 over Real Field with
1579 53 bits of precision
1583 The dimension of this algebra is `n^2`::
1585 sage: set_random_seed()
1586 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1587 sage: n = ZZ.random_element(1, n_max)
1588 sage: J = ComplexHermitianEJA(n)
1589 sage: J.dimension() == n^2
1592 The Jordan multiplication is what we think it is::
1594 sage: set_random_seed()
1595 sage: J = ComplexHermitianEJA.random_instance()
1596 sage: x,y = J.random_elements(2)
1597 sage: actual = (x*y).natural_representation()
1598 sage: X = x.natural_representation()
1599 sage: Y = y.natural_representation()
1600 sage: expected = (X*Y + Y*X)/2
1601 sage: actual == expected
1603 sage: J(expected) == x*y
1606 We can change the generator prefix::
1608 sage: ComplexHermitianEJA(2, prefix='z').gens()
1611 Our natural basis is normalized with respect to the natural inner
1612 product unless we specify otherwise::
1614 sage: set_random_seed()
1615 sage: J = ComplexHermitianEJA.random_instance()
1616 sage: all( b.norm() == 1 for b in J.gens() )
1619 Since our natural basis is normalized with respect to the natural
1620 inner product, and since we know that this algebra is an EJA, any
1621 left-multiplication operator's matrix will be symmetric because
1622 natural->EJA basis representation is an isometry and within the EJA
1623 the operator is self-adjoint by the Jordan axiom::
1625 sage: set_random_seed()
1626 sage: x = ComplexHermitianEJA.random_instance().random_element()
1627 sage: x.operator().matrix().is_symmetric()
1630 We can construct the (trivial) algebra of rank zero::
1632 sage: ComplexHermitianEJA(0)
1633 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1638 def _denormalized_basis(cls
, n
, field
):
1640 Returns a basis for the space of complex Hermitian n-by-n matrices.
1642 Why do we embed these? Basically, because all of numerical linear
1643 algebra assumes that you're working with vectors consisting of `n`
1644 entries from a field and scalars from the same field. There's no way
1645 to tell SageMath that (for example) the vectors contain complex
1646 numbers, while the scalar field is real.
1650 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1654 sage: set_random_seed()
1655 sage: n = ZZ.random_element(1,5)
1656 sage: field = QuadraticField(2, 'sqrt2')
1657 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1658 sage: all( M.is_symmetric() for M in B)
1662 R
= PolynomialRing(field
, 'z')
1664 F
= field
.extension(z
**2 + 1, 'I')
1667 # This is like the symmetric case, but we need to be careful:
1669 # * We want conjugate-symmetry, not just symmetry.
1670 # * The diagonal will (as a result) be real.
1674 for j
in range(i
+1):
1675 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1677 Sij
= cls
.real_embed(Eij
)
1680 # The second one has a minus because it's conjugated.
1681 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1683 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1686 # Since we embedded these, we can drop back to the "field" that we
1687 # started with instead of the complex extension "F".
1688 return ( s
.change_ring(field
) for s
in S
)
1691 def __init__(self
, n
, field
=AA
, **kwargs
):
1692 basis
= self
._denormalized
_basis
(n
,field
)
1693 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1694 self
.rank
.set_cache(n
)
1697 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1701 Embed the n-by-n quaternion matrix ``M`` into the space of real
1702 matrices of size 4n-by-4n by first sending each quaternion entry `z
1703 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1704 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1709 sage: from mjo.eja.eja_algebra import \
1710 ....: QuaternionMatrixEuclideanJordanAlgebra
1714 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1715 sage: i,j,k = Q.gens()
1716 sage: x = 1 + 2*i + 3*j + 4*k
1717 sage: M = matrix(Q, 1, [[x]])
1718 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1724 Embedding is a homomorphism (isomorphism, in fact)::
1726 sage: set_random_seed()
1727 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1728 sage: n = ZZ.random_element(n_max)
1729 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1730 sage: X = random_matrix(Q, n)
1731 sage: Y = random_matrix(Q, n)
1732 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1733 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1734 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1739 quaternions
= M
.base_ring()
1742 raise ValueError("the matrix 'M' must be square")
1744 F
= QuadraticField(-1, 'I')
1749 t
= z
.coefficient_tuple()
1754 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1755 [-c
+ d
*i
, a
- b
*i
]])
1756 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1757 blocks
.append(realM
)
1759 # We should have real entries by now, so use the realest field
1760 # we've got for the return value.
1761 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1766 def real_unembed(M
):
1768 The inverse of _embed_quaternion_matrix().
1772 sage: from mjo.eja.eja_algebra import \
1773 ....: QuaternionMatrixEuclideanJordanAlgebra
1777 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1778 ....: [-2, 1, -4, 3],
1779 ....: [-3, 4, 1, -2],
1780 ....: [-4, -3, 2, 1]])
1781 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1782 [1 + 2*i + 3*j + 4*k]
1786 Unembedding is the inverse of embedding::
1788 sage: set_random_seed()
1789 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1790 sage: M = random_matrix(Q, 3)
1791 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1792 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1798 raise ValueError("the matrix 'M' must be square")
1799 if not n
.mod(4).is_zero():
1800 raise ValueError("the matrix 'M' must be a quaternion embedding")
1802 # Use the base ring of the matrix to ensure that its entries can be
1803 # multiplied by elements of the quaternion algebra.
1804 field
= M
.base_ring()
1805 Q
= QuaternionAlgebra(field
,-1,-1)
1808 # Go top-left to bottom-right (reading order), converting every
1809 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1812 for l
in range(n
/4):
1813 for m
in range(n
/4):
1814 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1815 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1816 if submat
[0,0] != submat
[1,1].conjugate():
1817 raise ValueError('bad on-diagonal submatrix')
1818 if submat
[0,1] != -submat
[1,0].conjugate():
1819 raise ValueError('bad off-diagonal submatrix')
1820 z
= submat
[0,0].real()
1821 z
+= submat
[0,0].imag()*i
1822 z
+= submat
[0,1].real()*j
1823 z
+= submat
[0,1].imag()*k
1826 return matrix(Q
, n
/4, elements
)
1830 def natural_inner_product(cls
,X
,Y
):
1832 Compute a natural inner product in this algebra directly from
1837 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1841 This gives the same answer as the slow, default method implemented
1842 in :class:`MatrixEuclideanJordanAlgebra`::
1844 sage: set_random_seed()
1845 sage: J = QuaternionHermitianEJA.random_instance()
1846 sage: x,y = J.random_elements(2)
1847 sage: Xe = x.natural_representation()
1848 sage: Ye = y.natural_representation()
1849 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1850 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1851 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1852 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1853 sage: actual == expected
1857 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1860 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1862 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1863 matrices, the usual symmetric Jordan product, and the
1864 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1869 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1873 In theory, our "field" can be any subfield of the reals::
1875 sage: QuaternionHermitianEJA(2, RDF)
1876 Euclidean Jordan algebra of dimension 6 over Real Double Field
1877 sage: QuaternionHermitianEJA(2, RR)
1878 Euclidean Jordan algebra of dimension 6 over Real Field with
1879 53 bits of precision
1883 The dimension of this algebra is `2*n^2 - n`::
1885 sage: set_random_seed()
1886 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1887 sage: n = ZZ.random_element(1, n_max)
1888 sage: J = QuaternionHermitianEJA(n)
1889 sage: J.dimension() == 2*(n^2) - n
1892 The Jordan multiplication is what we think it is::
1894 sage: set_random_seed()
1895 sage: J = QuaternionHermitianEJA.random_instance()
1896 sage: x,y = J.random_elements(2)
1897 sage: actual = (x*y).natural_representation()
1898 sage: X = x.natural_representation()
1899 sage: Y = y.natural_representation()
1900 sage: expected = (X*Y + Y*X)/2
1901 sage: actual == expected
1903 sage: J(expected) == x*y
1906 We can change the generator prefix::
1908 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1909 (a0, a1, a2, a3, a4, a5)
1911 Our natural basis is normalized with respect to the natural inner
1912 product unless we specify otherwise::
1914 sage: set_random_seed()
1915 sage: J = QuaternionHermitianEJA.random_instance()
1916 sage: all( b.norm() == 1 for b in J.gens() )
1919 Since our natural basis is normalized with respect to the natural
1920 inner product, and since we know that this algebra is an EJA, any
1921 left-multiplication operator's matrix will be symmetric because
1922 natural->EJA basis representation is an isometry and within the EJA
1923 the operator is self-adjoint by the Jordan axiom::
1925 sage: set_random_seed()
1926 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1927 sage: x.operator().matrix().is_symmetric()
1930 We can construct the (trivial) algebra of rank zero::
1932 sage: QuaternionHermitianEJA(0)
1933 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1937 def _denormalized_basis(cls
, n
, field
):
1939 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1941 Why do we embed these? Basically, because all of numerical
1942 linear algebra assumes that you're working with vectors consisting
1943 of `n` entries from a field and scalars from the same field. There's
1944 no way to tell SageMath that (for example) the vectors contain
1945 complex numbers, while the scalar field is real.
1949 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1953 sage: set_random_seed()
1954 sage: n = ZZ.random_element(1,5)
1955 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1956 sage: all( M.is_symmetric() for M in B )
1960 Q
= QuaternionAlgebra(QQ
,-1,-1)
1963 # This is like the symmetric case, but we need to be careful:
1965 # * We want conjugate-symmetry, not just symmetry.
1966 # * The diagonal will (as a result) be real.
1970 for j
in range(i
+1):
1971 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1973 Sij
= cls
.real_embed(Eij
)
1976 # The second, third, and fourth ones have a minus
1977 # because they're conjugated.
1978 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1980 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1982 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1984 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1987 # Since we embedded these, we can drop back to the "field" that we
1988 # started with instead of the quaternion algebra "Q".
1989 return ( s
.change_ring(field
) for s
in S
)
1992 def __init__(self
, n
, field
=AA
, **kwargs
):
1993 basis
= self
._denormalized
_basis
(n
,field
)
1994 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1995 self
.rank
.set_cache(n
)
1998 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2000 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2001 with the half-trace inner product and jordan product ``x*y =
2002 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
2003 symmetric positive-definite "bilinear form" matrix. It has
2004 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
2005 when ``B`` is the identity matrix of order ``n-1``.
2009 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2010 ....: JordanSpinEJA)
2014 When no bilinear form is specified, the identity matrix is used,
2015 and the resulting algebra is the Jordan spin algebra::
2017 sage: J0 = BilinearFormEJA(3)
2018 sage: J1 = JordanSpinEJA(3)
2019 sage: J0.multiplication_table() == J0.multiplication_table()
2024 We can create a zero-dimensional algebra::
2026 sage: J = BilinearFormEJA(0)
2030 We can check the multiplication condition given in the Jordan, von
2031 Neumann, and Wigner paper (and also discussed on my "On the
2032 symmetry..." paper). Note that this relies heavily on the standard
2033 choice of basis, as does anything utilizing the bilinear form matrix::
2035 sage: set_random_seed()
2036 sage: n = ZZ.random_element(5)
2037 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2038 sage: B = M.transpose()*M
2039 sage: J = BilinearFormEJA(n, B=B)
2040 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2041 sage: V = J.vector_space()
2042 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2043 ....: for ei in eis ]
2044 sage: actual = [ sis[i]*sis[j]
2045 ....: for i in range(n-1)
2046 ....: for j in range(n-1) ]
2047 sage: expected = [ J.one() if i == j else J.zero()
2048 ....: for i in range(n-1)
2049 ....: for j in range(n-1) ]
2050 sage: actual == expected
2053 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
2055 self
._B
= matrix
.identity(field
, max(0,n
-1))
2059 V
= VectorSpace(field
, n
)
2060 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2069 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
2070 zbar
= y0
*xbar
+ x0
*ybar
2071 z
= V([z0
] + zbar
.list())
2072 mult_table
[i
][j
] = z
2074 # The rank of this algebra is two, unless we're in a
2075 # one-dimensional ambient space (because the rank is bounded
2076 # by the ambient dimension).
2077 fdeja
= super(BilinearFormEJA
, self
)
2078 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2079 self
.rank
.set_cache(min(n
,2))
2081 def inner_product(self
, x
, y
):
2083 Half of the trace inner product.
2085 This is defined so that the special case of the Jordan spin
2086 algebra gets the usual inner product.
2090 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2094 Ensure that this is one-half of the trace inner-product when
2095 the algebra isn't just the reals (when ``n`` isn't one). This
2096 is in Faraut and Koranyi, and also my "On the symmetry..."
2099 sage: set_random_seed()
2100 sage: n = ZZ.random_element(2,5)
2101 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2102 sage: B = M.transpose()*M
2103 sage: J = BilinearFormEJA(n, B=B)
2104 sage: x = J.random_element()
2105 sage: y = J.random_element()
2106 sage: x.inner_product(y) == (x*y).trace()/2
2110 xvec
= x
.to_vector()
2112 yvec
= y
.to_vector()
2114 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
2117 class JordanSpinEJA(BilinearFormEJA
):
2119 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2120 with the usual inner product and jordan product ``x*y =
2121 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2126 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2130 This multiplication table can be verified by hand::
2132 sage: J = JordanSpinEJA(4)
2133 sage: e0,e1,e2,e3 = J.gens()
2149 We can change the generator prefix::
2151 sage: JordanSpinEJA(2, prefix='B').gens()
2156 Ensure that we have the usual inner product on `R^n`::
2158 sage: set_random_seed()
2159 sage: J = JordanSpinEJA.random_instance()
2160 sage: x,y = J.random_elements(2)
2161 sage: X = x.natural_representation()
2162 sage: Y = y.natural_representation()
2163 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2167 def __init__(self
, n
, field
=AA
, **kwargs
):
2168 # This is a special case of the BilinearFormEJA with the identity
2169 # matrix as its bilinear form.
2170 return super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2173 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2175 The trivial Euclidean Jordan algebra consisting of only a zero element.
2179 sage: from mjo.eja.eja_algebra import TrivialEJA
2183 sage: J = TrivialEJA()
2190 sage: 7*J.one()*12*J.one()
2192 sage: J.one().inner_product(J.one())
2194 sage: J.one().norm()
2196 sage: J.one().subalgebra_generated_by()
2197 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2202 def __init__(self
, field
=AA
, **kwargs
):
2204 fdeja
= super(TrivialEJA
, self
)
2205 # The rank is zero using my definition, namely the dimension of the
2206 # largest subalgebra generated by any element.
2207 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2208 self
.rank
.set_cache(0)
2211 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2213 The external (orthogonal) direct sum of two other Euclidean Jordan
2214 algebras. Essentially the Cartesian product of its two factors.
2215 Every Euclidean Jordan algebra decomposes into an orthogonal
2216 direct sum of simple Euclidean Jordan algebras, so no generality
2217 is lost by providing only this construction.
2221 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2222 ....: RealSymmetricEJA,
2227 sage: J1 = HadamardEJA(2)
2228 sage: J2 = RealSymmetricEJA(3)
2229 sage: J = DirectSumEJA(J1,J2)
2236 def __init__(self
, J1
, J2
, field
=AA
, **kwargs
):
2240 V
= VectorSpace(field
, n
)
2241 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2245 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2246 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2250 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2251 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2253 fdeja
= super(DirectSumEJA
, self
)
2254 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2255 self
.rank
.set_cache(J1
.rank() + J2
.rank())