2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 from itertools
import repeat
10 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
11 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
12 from sage
.combinat
.free_module
import CombinatorialFreeModule
13 from sage
.matrix
.constructor
import matrix
14 from sage
.matrix
.matrix_space
import MatrixSpace
15 from sage
.misc
.cachefunc
import cached_method
16 from sage
.misc
.lazy_import
import lazy_import
17 from sage
.misc
.prandom
import choice
18 from sage
.misc
.table
import table
19 from sage
.modules
.free_module
import FreeModule
, VectorSpace
20 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
23 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
24 lazy_import('mjo.eja.eja_subalgebra',
25 'FiniteDimensionalEuclideanJordanSubalgebra')
26 from mjo
.eja
.eja_utils
import _mat2vec
28 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
30 def _coerce_map_from_base_ring(self
):
32 Disable the map from the base ring into the algebra.
34 Performing a nonsense conversion like this automatically
35 is counterpedagogical. The fallback is to try the usual
36 element constructor, which should also fail.
40 sage: from mjo.eja.eja_algebra import random_eja
44 sage: set_random_seed()
45 sage: J = random_eja()
47 Traceback (most recent call last):
49 ValueError: not a naturally-represented algebra element
64 sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
68 By definition, Jordan multiplication commutes::
70 sage: set_random_seed()
71 sage: J = random_eja()
72 sage: x,y = J.random_elements(2)
78 The ``field`` we're given must be real::
80 sage: JordanSpinEJA(2,QQbar)
81 Traceback (most recent call last):
83 ValueError: field is not real
87 if not field
.is_subring(RR
):
88 # Note: this does return true for the real algebraic
89 # field, and any quadratic field where we've specified
91 raise ValueError('field is not real')
93 self
._natural
_basis
= natural_basis
96 category
= MagmaticAlgebras(field
).FiniteDimensional()
97 category
= category
.WithBasis().Unital()
99 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
101 range(len(mult_table
)),
104 self
.print_options(bracket
='')
106 # The multiplication table we're given is necessarily in terms
107 # of vectors, because we don't have an algebra yet for
108 # anything to be an element of. However, it's faster in the
109 # long run to have the multiplication table be in terms of
110 # algebra elements. We do this after calling the superclass
111 # constructor so that from_vector() knows what to do.
112 self
._multiplication
_table
= [
113 list(map(lambda x
: self
.from_vector(x
), ls
))
118 def _element_constructor_(self
, elt
):
120 Construct an element of this algebra from its natural
123 This gets called only after the parent element _call_ method
124 fails to find a coercion for the argument.
128 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
130 ....: RealSymmetricEJA)
134 The identity in `S^n` is converted to the identity in the EJA::
136 sage: J = RealSymmetricEJA(3)
137 sage: I = matrix.identity(QQ,3)
138 sage: J(I) == J.one()
141 This skew-symmetric matrix can't be represented in the EJA::
143 sage: J = RealSymmetricEJA(3)
144 sage: A = matrix(QQ,3, lambda i,j: i-j)
146 Traceback (most recent call last):
148 ArithmeticError: vector is not in free module
152 Ensure that we can convert any element of the two non-matrix
153 simple algebras (whose natural representations are their usual
154 vector representations) back and forth faithfully::
156 sage: set_random_seed()
157 sage: J = HadamardEJA.random_instance()
158 sage: x = J.random_element()
159 sage: J(x.to_vector().column()) == x
161 sage: J = JordanSpinEJA.random_instance()
162 sage: x = J.random_element()
163 sage: J(x.to_vector().column()) == x
167 msg
= "not a naturally-represented algebra element"
169 # The superclass implementation of random_element()
170 # needs to be able to coerce "0" into the algebra.
172 elif elt
in self
.base_ring():
173 # Ensure that no base ring -> algebra coercion is performed
174 # by this method. There's some stupidity in sage that would
175 # otherwise propagate to this method; for example, sage thinks
176 # that the integer 3 belongs to the space of 2-by-2 matrices.
177 raise ValueError(msg
)
179 natural_basis
= self
.natural_basis()
180 basis_space
= natural_basis
[0].matrix_space()
181 if elt
not in basis_space
:
182 raise ValueError(msg
)
184 # Thanks for nothing! Matrix spaces aren't vector spaces in
185 # Sage, so we have to figure out its natural-basis coordinates
186 # ourselves. We use the basis space's ring instead of the
187 # element's ring because the basis space might be an algebraic
188 # closure whereas the base ring of the 3-by-3 identity matrix
189 # could be QQ instead of QQbar.
190 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
191 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
192 coords
= W
.coordinate_vector(_mat2vec(elt
))
193 return self
.from_vector(coords
)
196 def _max_test_case_size():
198 Return an integer "size" that is an upper bound on the size of
199 this algebra when it is used in a random test
200 case. Unfortunately, the term "size" is quite vague -- when
201 dealing with `R^n` under either the Hadamard or Jordan spin
202 product, the "size" refers to the dimension `n`. When dealing
203 with a matrix algebra (real symmetric or complex/quaternion
204 Hermitian), it refers to the size of the matrix, which is
205 far less than the dimension of the underlying vector space.
207 We default to five in this class, which is safe in `R^n`. The
208 matrix algebra subclasses (or any class where the "size" is
209 interpreted to be far less than the dimension) should override
210 with a smaller number.
216 Return a string representation of ``self``.
220 sage: from mjo.eja.eja_algebra import JordanSpinEJA
224 Ensure that it says what we think it says::
226 sage: JordanSpinEJA(2, field=AA)
227 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
228 sage: JordanSpinEJA(3, field=RDF)
229 Euclidean Jordan algebra of dimension 3 over Real Double Field
232 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
233 return fmt
.format(self
.dimension(), self
.base_ring())
235 def product_on_basis(self
, i
, j
):
236 return self
._multiplication
_table
[i
][j
]
239 def characteristic_polynomial_of(self
):
241 Return the algebra's "characteristic polynomial of" function,
242 which is itself a multivariate polynomial that, when evaluated
243 at the coordinates of some algebra element, returns that
244 element's characteristic polynomial.
246 The resulting polynomial has `n+1` variables, where `n` is the
247 dimension of this algebra. The first `n` variables correspond to
248 the coordinates of an algebra element: when evaluated at the
249 coordinates of an algebra element with respect to a certain
250 basis, the result is a univariate polynomial (in the one
251 remaining variable ``t``), namely the characteristic polynomial
256 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
260 The characteristic polynomial in the spin algebra is given in
261 Alizadeh, Example 11.11::
263 sage: J = JordanSpinEJA(3)
264 sage: p = J.characteristic_polynomial_of(); p
265 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
266 sage: xvec = J.one().to_vector()
270 By definition, the characteristic polynomial is a monic
271 degree-zero polynomial in a rank-zero algebra. Note that
272 Cayley-Hamilton is indeed satisfied since the polynomial
273 ``1`` evaluates to the identity element of the algebra on
276 sage: J = TrivialEJA()
277 sage: J.characteristic_polynomial_of()
284 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
285 a
= self
._charpoly
_coefficients
()
287 # We go to a bit of trouble here to reorder the
288 # indeterminates, so that it's easier to evaluate the
289 # characteristic polynomial at x's coordinates and get back
290 # something in terms of t, which is what we want.
291 S
= PolynomialRing(self
.base_ring(),'t')
295 S
= PolynomialRing(S
, R
.variable_names())
298 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
301 def inner_product(self
, x
, y
):
303 The inner product associated with this Euclidean Jordan algebra.
305 Defaults to the trace inner product, but can be overridden by
306 subclasses if they are sure that the necessary properties are
311 sage: from mjo.eja.eja_algebra import random_eja
315 Our inner product is "associative," which means the following for
316 a symmetric bilinear form::
318 sage: set_random_seed()
319 sage: J = random_eja()
320 sage: x,y,z = J.random_elements(3)
321 sage: (x*y).inner_product(z) == y.inner_product(x*z)
325 X
= x
.natural_representation()
326 Y
= y
.natural_representation()
327 return self
.natural_inner_product(X
,Y
)
330 def is_trivial(self
):
332 Return whether or not this algebra is trivial.
334 A trivial algebra contains only the zero element.
338 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
343 sage: J = ComplexHermitianEJA(3)
349 sage: J = TrivialEJA()
354 return self
.dimension() == 0
357 def multiplication_table(self
):
359 Return a visual representation of this algebra's multiplication
360 table (on basis elements).
364 sage: from mjo.eja.eja_algebra import JordanSpinEJA
368 sage: J = JordanSpinEJA(4)
369 sage: J.multiplication_table()
370 +----++----+----+----+----+
371 | * || e0 | e1 | e2 | e3 |
372 +====++====+====+====+====+
373 | e0 || e0 | e1 | e2 | e3 |
374 +----++----+----+----+----+
375 | e1 || e1 | e0 | 0 | 0 |
376 +----++----+----+----+----+
377 | e2 || e2 | 0 | e0 | 0 |
378 +----++----+----+----+----+
379 | e3 || e3 | 0 | 0 | e0 |
380 +----++----+----+----+----+
383 M
= list(self
._multiplication
_table
) # copy
384 for i
in range(len(M
)):
385 # M had better be "square"
386 M
[i
] = [self
.monomial(i
)] + M
[i
]
387 M
= [["*"] + list(self
.gens())] + M
388 return table(M
, header_row
=True, header_column
=True, frame
=True)
391 def natural_basis(self
):
393 Return a more-natural representation of this algebra's basis.
395 Every finite-dimensional Euclidean Jordan Algebra is a direct
396 sum of five simple algebras, four of which comprise Hermitian
397 matrices. This method returns the original "natural" basis
398 for our underlying vector space. (Typically, the natural basis
399 is used to construct the multiplication table in the first place.)
401 Note that this will always return a matrix. The standard basis
402 in `R^n` will be returned as `n`-by-`1` column matrices.
406 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
407 ....: RealSymmetricEJA)
411 sage: J = RealSymmetricEJA(2)
413 Finite family {0: e0, 1: e1, 2: e2}
414 sage: J.natural_basis()
416 [1 0] [ 0 0.7071067811865475?] [0 0]
417 [0 0], [0.7071067811865475? 0], [0 1]
422 sage: J = JordanSpinEJA(2)
424 Finite family {0: e0, 1: e1}
425 sage: J.natural_basis()
432 if self
._natural
_basis
is None:
433 M
= self
.natural_basis_space()
434 return tuple( M(b
.to_vector()) for b
in self
.basis() )
436 return self
._natural
_basis
439 def natural_basis_space(self
):
441 Return the matrix space in which this algebra's natural basis
444 Generally this will be an `n`-by-`1` column-vector space,
445 except when the algebra is trivial. There it's `n`-by-`n`
446 (where `n` is zero), to ensure that two elements of the
447 natural basis space (empty matrices) can be multiplied.
449 if self
.is_trivial():
450 return MatrixSpace(self
.base_ring(), 0)
451 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
452 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
454 return self
._natural
_basis
[0].matrix_space()
458 def natural_inner_product(X
,Y
):
460 Compute the inner product of two naturally-represented elements.
462 For example in the real symmetric matrix EJA, this will compute
463 the trace inner-product of two n-by-n symmetric matrices. The
464 default should work for the real cartesian product EJA, the
465 Jordan spin EJA, and the real symmetric matrices. The others
466 will have to be overridden.
468 return (X
.conjugate_transpose()*Y
).trace()
474 Return the unit element of this algebra.
478 sage: from mjo.eja.eja_algebra import (HadamardEJA,
483 sage: J = HadamardEJA(5)
485 e0 + e1 + e2 + e3 + e4
489 The identity element acts like the identity::
491 sage: set_random_seed()
492 sage: J = random_eja()
493 sage: x = J.random_element()
494 sage: J.one()*x == x and x*J.one() == x
497 The matrix of the unit element's operator is the identity::
499 sage: set_random_seed()
500 sage: J = random_eja()
501 sage: actual = J.one().operator().matrix()
502 sage: expected = matrix.identity(J.base_ring(), J.dimension())
503 sage: actual == expected
507 # We can brute-force compute the matrices of the operators
508 # that correspond to the basis elements of this algebra.
509 # If some linear combination of those basis elements is the
510 # algebra identity, then the same linear combination of
511 # their matrices has to be the identity matrix.
513 # Of course, matrices aren't vectors in sage, so we have to
514 # appeal to the "long vectors" isometry.
515 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
517 # Now we use basis linear algebra to find the coefficients,
518 # of the matrices-as-vectors-linear-combination, which should
519 # work for the original algebra basis too.
520 A
= matrix
.column(self
.base_ring(), oper_vecs
)
522 # We used the isometry on the left-hand side already, but we
523 # still need to do it for the right-hand side. Recall that we
524 # wanted something that summed to the identity matrix.
525 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
527 # Now if there's an identity element in the algebra, this should work.
528 coeffs
= A
.solve_right(b
)
529 return self
.linear_combination(zip(self
.gens(), coeffs
))
532 def peirce_decomposition(self
, c
):
534 The Peirce decomposition of this algebra relative to the
537 In the future, this can be extended to a complete system of
538 orthogonal idempotents.
542 - ``c`` -- an idempotent of this algebra.
546 A triple (J0, J5, J1) containing two subalgebras and one subspace
549 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
550 corresponding to the eigenvalue zero.
552 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
553 corresponding to the eigenvalue one-half.
555 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
556 corresponding to the eigenvalue one.
558 These are the only possible eigenspaces for that operator, and this
559 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
560 orthogonal, and are subalgebras of this algebra with the appropriate
565 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
569 The canonical example comes from the symmetric matrices, which
570 decompose into diagonal and off-diagonal parts::
572 sage: J = RealSymmetricEJA(3)
573 sage: C = matrix(QQ, [ [1,0,0],
577 sage: J0,J5,J1 = J.peirce_decomposition(c)
579 Euclidean Jordan algebra of dimension 1...
581 Vector space of degree 6 and dimension 2...
583 Euclidean Jordan algebra of dimension 3...
587 Every algebra decomposes trivially with respect to its identity
590 sage: set_random_seed()
591 sage: J = random_eja()
592 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
593 sage: J0.dimension() == 0 and J5.dimension() == 0
595 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
598 The identity elements in the two subalgebras are the
599 projections onto their respective subspaces of the
600 superalgebra's identity element::
602 sage: set_random_seed()
603 sage: J = random_eja()
604 sage: x = J.random_element()
605 sage: if not J.is_trivial():
606 ....: while x.is_nilpotent():
607 ....: x = J.random_element()
608 sage: c = x.subalgebra_idempotent()
609 sage: J0,J5,J1 = J.peirce_decomposition(c)
610 sage: J1(c) == J1.one()
612 sage: J0(J.one() - c) == J0.one()
616 if not c
.is_idempotent():
617 raise ValueError("element is not idempotent: %s" % c
)
619 # Default these to what they should be if they turn out to be
620 # trivial, because eigenspaces_left() won't return eigenvalues
621 # corresponding to trivial spaces (e.g. it returns only the
622 # eigenspace corresponding to lambda=1 if you take the
623 # decomposition relative to the identity element).
624 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
625 J0
= trivial
# eigenvalue zero
626 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
627 J1
= trivial
# eigenvalue one
629 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
630 if eigval
== ~
(self
.base_ring()(2)):
633 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
634 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
640 raise ValueError("unexpected eigenvalue: %s" % eigval
)
645 def random_elements(self
, count
):
647 Return ``count`` random elements as a tuple.
651 sage: from mjo.eja.eja_algebra import JordanSpinEJA
655 sage: J = JordanSpinEJA(3)
656 sage: x,y,z = J.random_elements(3)
657 sage: all( [ x in J, y in J, z in J ])
659 sage: len( J.random_elements(10) ) == 10
663 return tuple( self
.random_element() for idx
in range(count
) )
666 def random_instance(cls
, field
=AA
, **kwargs
):
668 Return a random instance of this type of algebra.
670 Beware, this will crash for "most instances" because the
671 constructor below looks wrong.
673 if cls
is TrivialEJA
:
674 # The TrivialEJA class doesn't take an "n" argument because
678 n
= ZZ
.random_element(cls
._max
_test
_case
_size
() + 1)
679 return cls(n
, field
, **kwargs
)
682 def _charpoly_coefficients(self
):
684 The `r` polynomial coefficients of the "characteristic polynomial
688 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
689 R
= PolynomialRing(self
.base_ring(), var_names
)
691 F
= R
.fraction_field()
694 # From a result in my book, these are the entries of the
695 # basis representation of L_x.
696 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
699 L_x
= matrix(F
, n
, n
, L_x_i_j
)
702 if self
.rank
.is_in_cache():
704 # There's no need to pad the system with redundant
705 # columns if we *know* they'll be redundant.
708 # Compute an extra power in case the rank is equal to
709 # the dimension (otherwise, we would stop at x^(r-1)).
710 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
711 for k
in range(n
+1) ]
712 A
= matrix
.column(F
, x_powers
[:n
])
713 AE
= A
.extended_echelon_form()
720 # The theory says that only the first "r" coefficients are
721 # nonzero, and they actually live in the original polynomial
722 # ring and not the fraction field. We negate them because
723 # in the actual characteristic polynomial, they get moved
724 # to the other side where x^r lives.
725 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
730 Return the rank of this EJA.
732 This is a cached method because we know the rank a priori for
733 all of the algebras we can construct. Thus we can avoid the
734 expensive ``_charpoly_coefficients()`` call unless we truly
735 need to compute the whole characteristic polynomial.
739 sage: from mjo.eja.eja_algebra import (HadamardEJA,
741 ....: RealSymmetricEJA,
742 ....: ComplexHermitianEJA,
743 ....: QuaternionHermitianEJA,
748 The rank of the Jordan spin algebra is always two::
750 sage: JordanSpinEJA(2).rank()
752 sage: JordanSpinEJA(3).rank()
754 sage: JordanSpinEJA(4).rank()
757 The rank of the `n`-by-`n` Hermitian real, complex, or
758 quaternion matrices is `n`::
760 sage: RealSymmetricEJA(4).rank()
762 sage: ComplexHermitianEJA(3).rank()
764 sage: QuaternionHermitianEJA(2).rank()
769 Ensure that every EJA that we know how to construct has a
770 positive integer rank, unless the algebra is trivial in
771 which case its rank will be zero::
773 sage: set_random_seed()
774 sage: J = random_eja()
778 sage: r > 0 or (r == 0 and J.is_trivial())
781 Ensure that computing the rank actually works, since the ranks
782 of all simple algebras are known and will be cached by default::
784 sage: J = HadamardEJA(4)
785 sage: J.rank.clear_cache()
791 sage: J = JordanSpinEJA(4)
792 sage: J.rank.clear_cache()
798 sage: J = RealSymmetricEJA(3)
799 sage: J.rank.clear_cache()
805 sage: J = ComplexHermitianEJA(2)
806 sage: J.rank.clear_cache()
812 sage: J = QuaternionHermitianEJA(2)
813 sage: J.rank.clear_cache()
817 return len(self
._charpoly
_coefficients
())
820 def vector_space(self
):
822 Return the vector space that underlies this algebra.
826 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
830 sage: J = RealSymmetricEJA(2)
831 sage: J.vector_space()
832 Vector space of dimension 3 over...
835 return self
.zero().to_vector().parent().ambient_vector_space()
838 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
841 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
):
843 Return the Euclidean Jordan Algebra corresponding to the set
844 `R^n` under the Hadamard product.
846 Note: this is nothing more than the Cartesian product of ``n``
847 copies of the spin algebra. Once Cartesian product algebras
848 are implemented, this can go.
852 sage: from mjo.eja.eja_algebra import HadamardEJA
856 This multiplication table can be verified by hand::
858 sage: J = HadamardEJA(3)
859 sage: e0,e1,e2 = J.gens()
875 We can change the generator prefix::
877 sage: HadamardEJA(3, prefix='r').gens()
881 def __init__(self
, n
, field
=AA
, **kwargs
):
882 V
= VectorSpace(field
, n
)
883 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
886 fdeja
= super(HadamardEJA
, self
)
887 fdeja
.__init
__(field
, mult_table
, **kwargs
)
888 self
.rank
.set_cache(n
)
890 def inner_product(self
, x
, y
):
892 Faster to reimplement than to use natural representations.
896 sage: from mjo.eja.eja_algebra import HadamardEJA
900 Ensure that this is the usual inner product for the algebras
903 sage: set_random_seed()
904 sage: J = HadamardEJA.random_instance()
905 sage: x,y = J.random_elements(2)
906 sage: X = x.natural_representation()
907 sage: Y = y.natural_representation()
908 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
912 return x
.to_vector().inner_product(y
.to_vector())
915 def random_eja(field
=AA
):
917 Return a "random" finite-dimensional Euclidean Jordan Algebra.
921 sage: from mjo.eja.eja_algebra import random_eja
926 Euclidean Jordan algebra of dimension...
929 classname
= choice([TrivialEJA
,
934 QuaternionHermitianEJA
])
935 return classname
.random_instance(field
=field
)
940 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
942 def _max_test_case_size():
943 # Play it safe, since this will be squared and the underlying
944 # field can have dimension 4 (quaternions) too.
947 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
949 Compared to the superclass constructor, we take a basis instead of
950 a multiplication table because the latter can be computed in terms
951 of the former when the product is known (like it is here).
953 # Used in this class's fast _charpoly_coefficients() override.
954 self
._basis
_normalizers
= None
956 # We're going to loop through this a few times, so now's a good
957 # time to ensure that it isn't a generator expression.
960 if len(basis
) > 1 and normalize_basis
:
961 # We'll need sqrt(2) to normalize the basis, and this
962 # winds up in the multiplication table, so the whole
963 # algebra needs to be over the field extension.
964 R
= PolynomialRing(field
, 'z')
967 if p
.is_irreducible():
968 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
969 basis
= tuple( s
.change_ring(field
) for s
in basis
)
970 self
._basis
_normalizers
= tuple(
971 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
972 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
974 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
976 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
977 fdeja
.__init
__(field
, Qs
, natural_basis
=basis
, **kwargs
)
982 def _charpoly_coefficients(self
):
984 Override the parent method with something that tries to compute
985 over a faster (non-extension) field.
987 if self
._basis
_normalizers
is None:
988 # We didn't normalize, so assume that the basis we started
989 # with had entries in a nice field.
990 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
992 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
993 self
._basis
_normalizers
) )
995 # Do this over the rationals and convert back at the end.
996 # Only works because we know the entries of the basis are
998 J
= MatrixEuclideanJordanAlgebra(QQ
,
1000 normalize_basis
=False)
1001 a
= J
._charpoly
_coefficients
()
1003 # Unfortunately, changing the basis does change the
1004 # coefficients of the characteristic polynomial, but since
1005 # these are really the coefficients of the "characteristic
1006 # polynomial of" function, everything is still nice and
1007 # unevaluated. It's therefore "obvious" how scaling the
1008 # basis affects the coordinate variables X1, X2, et
1009 # cetera. Scaling the first basis vector up by "n" adds a
1010 # factor of 1/n into every "X1" term, for example. So here
1011 # we simply undo the basis_normalizer scaling that we
1012 # performed earlier.
1014 # The a[0] access here is safe because trivial algebras
1015 # won't have any basis normalizers and therefore won't
1016 # make it to this "else" branch.
1017 XS
= a
[0].parent().gens()
1018 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1019 for i
in range(len(XS
)) }
1020 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1024 def multiplication_table_from_matrix_basis(basis
):
1026 At least three of the five simple Euclidean Jordan algebras have the
1027 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1028 multiplication on the right is matrix multiplication. Given a basis
1029 for the underlying matrix space, this function returns a
1030 multiplication table (obtained by looping through the basis
1031 elements) for an algebra of those matrices.
1033 # In S^2, for example, we nominally have four coordinates even
1034 # though the space is of dimension three only. The vector space V
1035 # is supposed to hold the entire long vector, and the subspace W
1036 # of V will be spanned by the vectors that arise from symmetric
1037 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1041 field
= basis
[0].base_ring()
1042 dimension
= basis
[0].nrows()
1044 V
= VectorSpace(field
, dimension
**2)
1045 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1047 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1050 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1051 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1059 Embed the matrix ``M`` into a space of real matrices.
1061 The matrix ``M`` can have entries in any field at the moment:
1062 the real numbers, complex numbers, or quaternions. And although
1063 they are not a field, we can probably support octonions at some
1064 point, too. This function returns a real matrix that "acts like"
1065 the original with respect to matrix multiplication; i.e.
1067 real_embed(M*N) = real_embed(M)*real_embed(N)
1070 raise NotImplementedError
1074 def real_unembed(M
):
1076 The inverse of :meth:`real_embed`.
1078 raise NotImplementedError
1082 def natural_inner_product(cls
,X
,Y
):
1083 Xu
= cls
.real_unembed(X
)
1084 Yu
= cls
.real_unembed(Y
)
1085 tr
= (Xu
*Yu
).trace()
1088 # It's real already.
1091 # Otherwise, try the thing that works for complex numbers; and
1092 # if that doesn't work, the thing that works for quaternions.
1094 return tr
.vector()[0] # real part, imag part is index 1
1095 except AttributeError:
1096 # A quaternions doesn't have a vector() method, but does
1097 # have coefficient_tuple() method that returns the
1098 # coefficients of 1, i, j, and k -- in that order.
1099 return tr
.coefficient_tuple()[0]
1102 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1106 The identity function, for embedding real matrices into real
1112 def real_unembed(M
):
1114 The identity function, for unembedding real matrices from real
1120 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1122 The rank-n simple EJA consisting of real symmetric n-by-n
1123 matrices, the usual symmetric Jordan product, and the trace inner
1124 product. It has dimension `(n^2 + n)/2` over the reals.
1128 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1132 sage: J = RealSymmetricEJA(2)
1133 sage: e0, e1, e2 = J.gens()
1141 In theory, our "field" can be any subfield of the reals::
1143 sage: RealSymmetricEJA(2, RDF)
1144 Euclidean Jordan algebra of dimension 3 over Real Double Field
1145 sage: RealSymmetricEJA(2, RR)
1146 Euclidean Jordan algebra of dimension 3 over Real Field with
1147 53 bits of precision
1151 The dimension of this algebra is `(n^2 + n) / 2`::
1153 sage: set_random_seed()
1154 sage: n_max = RealSymmetricEJA._max_test_case_size()
1155 sage: n = ZZ.random_element(1, n_max)
1156 sage: J = RealSymmetricEJA(n)
1157 sage: J.dimension() == (n^2 + n)/2
1160 The Jordan multiplication is what we think it is::
1162 sage: set_random_seed()
1163 sage: J = RealSymmetricEJA.random_instance()
1164 sage: x,y = J.random_elements(2)
1165 sage: actual = (x*y).natural_representation()
1166 sage: X = x.natural_representation()
1167 sage: Y = y.natural_representation()
1168 sage: expected = (X*Y + Y*X)/2
1169 sage: actual == expected
1171 sage: J(expected) == x*y
1174 We can change the generator prefix::
1176 sage: RealSymmetricEJA(3, prefix='q').gens()
1177 (q0, q1, q2, q3, q4, q5)
1179 Our natural basis is normalized with respect to the natural inner
1180 product unless we specify otherwise::
1182 sage: set_random_seed()
1183 sage: J = RealSymmetricEJA.random_instance()
1184 sage: all( b.norm() == 1 for b in J.gens() )
1187 Since our natural basis is normalized with respect to the natural
1188 inner product, and since we know that this algebra is an EJA, any
1189 left-multiplication operator's matrix will be symmetric because
1190 natural->EJA basis representation is an isometry and within the EJA
1191 the operator is self-adjoint by the Jordan axiom::
1193 sage: set_random_seed()
1194 sage: x = RealSymmetricEJA.random_instance().random_element()
1195 sage: x.operator().matrix().is_symmetric()
1198 We can construct the (trivial) algebra of rank zero::
1200 sage: RealSymmetricEJA(0)
1201 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1205 def _denormalized_basis(cls
, n
, field
):
1207 Return a basis for the space of real symmetric n-by-n matrices.
1211 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1215 sage: set_random_seed()
1216 sage: n = ZZ.random_element(1,5)
1217 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1218 sage: all( M.is_symmetric() for M in B)
1222 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1226 for j
in range(i
+1):
1227 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1231 Sij
= Eij
+ Eij
.transpose()
1237 def _max_test_case_size():
1238 return 4 # Dimension 10
1241 def __init__(self
, n
, field
=AA
, **kwargs
):
1242 basis
= self
._denormalized
_basis
(n
, field
)
1243 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, **kwargs
)
1244 self
.rank
.set_cache(n
)
1247 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1251 Embed the n-by-n complex matrix ``M`` into the space of real
1252 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1253 bi` to the block matrix ``[[a,b],[-b,a]]``.
1257 sage: from mjo.eja.eja_algebra import \
1258 ....: ComplexMatrixEuclideanJordanAlgebra
1262 sage: F = QuadraticField(-1, 'I')
1263 sage: x1 = F(4 - 2*i)
1264 sage: x2 = F(1 + 2*i)
1267 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1268 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1277 Embedding is a homomorphism (isomorphism, in fact)::
1279 sage: set_random_seed()
1280 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1281 sage: n = ZZ.random_element(n_max)
1282 sage: F = QuadraticField(-1, 'I')
1283 sage: X = random_matrix(F, n)
1284 sage: Y = random_matrix(F, n)
1285 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1286 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1287 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1294 raise ValueError("the matrix 'M' must be square")
1296 # We don't need any adjoined elements...
1297 field
= M
.base_ring().base_ring()
1301 a
= z
.list()[0] # real part, I guess
1302 b
= z
.list()[1] # imag part, I guess
1303 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1305 return matrix
.block(field
, n
, blocks
)
1309 def real_unembed(M
):
1311 The inverse of _embed_complex_matrix().
1315 sage: from mjo.eja.eja_algebra import \
1316 ....: ComplexMatrixEuclideanJordanAlgebra
1320 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1321 ....: [-2, 1, -4, 3],
1322 ....: [ 9, 10, 11, 12],
1323 ....: [-10, 9, -12, 11] ])
1324 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1326 [ 10*I + 9 12*I + 11]
1330 Unembedding is the inverse of embedding::
1332 sage: set_random_seed()
1333 sage: F = QuadraticField(-1, 'I')
1334 sage: M = random_matrix(F, 3)
1335 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1336 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1342 raise ValueError("the matrix 'M' must be square")
1343 if not n
.mod(2).is_zero():
1344 raise ValueError("the matrix 'M' must be a complex embedding")
1346 # If "M" was normalized, its base ring might have roots
1347 # adjoined and they can stick around after unembedding.
1348 field
= M
.base_ring()
1349 R
= PolynomialRing(field
, 'z')
1352 # Sage doesn't know how to embed AA into QQbar, i.e. how
1353 # to adjoin sqrt(-1) to AA.
1356 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1359 # Go top-left to bottom-right (reading order), converting every
1360 # 2-by-2 block we see to a single complex element.
1362 for k
in range(n
/2):
1363 for j
in range(n
/2):
1364 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1365 if submat
[0,0] != submat
[1,1]:
1366 raise ValueError('bad on-diagonal submatrix')
1367 if submat
[0,1] != -submat
[1,0]:
1368 raise ValueError('bad off-diagonal submatrix')
1369 z
= submat
[0,0] + submat
[0,1]*i
1372 return matrix(F
, n
/2, elements
)
1376 def natural_inner_product(cls
,X
,Y
):
1378 Compute a natural inner product in this algebra directly from
1383 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1387 This gives the same answer as the slow, default method implemented
1388 in :class:`MatrixEuclideanJordanAlgebra`::
1390 sage: set_random_seed()
1391 sage: J = ComplexHermitianEJA.random_instance()
1392 sage: x,y = J.random_elements(2)
1393 sage: Xe = x.natural_representation()
1394 sage: Ye = y.natural_representation()
1395 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1396 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1397 sage: expected = (X*Y).trace().real()
1398 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1399 sage: actual == expected
1403 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1406 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1408 The rank-n simple EJA consisting of complex Hermitian n-by-n
1409 matrices over the real numbers, the usual symmetric Jordan product,
1410 and the real-part-of-trace inner product. It has dimension `n^2` over
1415 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1419 In theory, our "field" can be any subfield of the reals::
1421 sage: ComplexHermitianEJA(2, RDF)
1422 Euclidean Jordan algebra of dimension 4 over Real Double Field
1423 sage: ComplexHermitianEJA(2, RR)
1424 Euclidean Jordan algebra of dimension 4 over Real Field with
1425 53 bits of precision
1429 The dimension of this algebra is `n^2`::
1431 sage: set_random_seed()
1432 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1433 sage: n = ZZ.random_element(1, n_max)
1434 sage: J = ComplexHermitianEJA(n)
1435 sage: J.dimension() == n^2
1438 The Jordan multiplication is what we think it is::
1440 sage: set_random_seed()
1441 sage: J = ComplexHermitianEJA.random_instance()
1442 sage: x,y = J.random_elements(2)
1443 sage: actual = (x*y).natural_representation()
1444 sage: X = x.natural_representation()
1445 sage: Y = y.natural_representation()
1446 sage: expected = (X*Y + Y*X)/2
1447 sage: actual == expected
1449 sage: J(expected) == x*y
1452 We can change the generator prefix::
1454 sage: ComplexHermitianEJA(2, prefix='z').gens()
1457 Our natural basis is normalized with respect to the natural inner
1458 product unless we specify otherwise::
1460 sage: set_random_seed()
1461 sage: J = ComplexHermitianEJA.random_instance()
1462 sage: all( b.norm() == 1 for b in J.gens() )
1465 Since our natural basis is normalized with respect to the natural
1466 inner product, and since we know that this algebra is an EJA, any
1467 left-multiplication operator's matrix will be symmetric because
1468 natural->EJA basis representation is an isometry and within the EJA
1469 the operator is self-adjoint by the Jordan axiom::
1471 sage: set_random_seed()
1472 sage: x = ComplexHermitianEJA.random_instance().random_element()
1473 sage: x.operator().matrix().is_symmetric()
1476 We can construct the (trivial) algebra of rank zero::
1478 sage: ComplexHermitianEJA(0)
1479 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1484 def _denormalized_basis(cls
, n
, field
):
1486 Returns a basis for the space of complex Hermitian n-by-n matrices.
1488 Why do we embed these? Basically, because all of numerical linear
1489 algebra assumes that you're working with vectors consisting of `n`
1490 entries from a field and scalars from the same field. There's no way
1491 to tell SageMath that (for example) the vectors contain complex
1492 numbers, while the scalar field is real.
1496 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1500 sage: set_random_seed()
1501 sage: n = ZZ.random_element(1,5)
1502 sage: field = QuadraticField(2, 'sqrt2')
1503 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1504 sage: all( M.is_symmetric() for M in B)
1508 R
= PolynomialRing(field
, 'z')
1510 F
= field
.extension(z
**2 + 1, 'I')
1513 # This is like the symmetric case, but we need to be careful:
1515 # * We want conjugate-symmetry, not just symmetry.
1516 # * The diagonal will (as a result) be real.
1520 for j
in range(i
+1):
1521 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1523 Sij
= cls
.real_embed(Eij
)
1526 # The second one has a minus because it's conjugated.
1527 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1529 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1532 # Since we embedded these, we can drop back to the "field" that we
1533 # started with instead of the complex extension "F".
1534 return ( s
.change_ring(field
) for s
in S
)
1537 def __init__(self
, n
, field
=AA
, **kwargs
):
1538 basis
= self
._denormalized
_basis
(n
,field
)
1539 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1540 self
.rank
.set_cache(n
)
1543 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1547 Embed the n-by-n quaternion matrix ``M`` into the space of real
1548 matrices of size 4n-by-4n by first sending each quaternion entry `z
1549 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1550 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1555 sage: from mjo.eja.eja_algebra import \
1556 ....: QuaternionMatrixEuclideanJordanAlgebra
1560 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1561 sage: i,j,k = Q.gens()
1562 sage: x = 1 + 2*i + 3*j + 4*k
1563 sage: M = matrix(Q, 1, [[x]])
1564 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1570 Embedding is a homomorphism (isomorphism, in fact)::
1572 sage: set_random_seed()
1573 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1574 sage: n = ZZ.random_element(n_max)
1575 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1576 sage: X = random_matrix(Q, n)
1577 sage: Y = random_matrix(Q, n)
1578 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1579 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1580 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1585 quaternions
= M
.base_ring()
1588 raise ValueError("the matrix 'M' must be square")
1590 F
= QuadraticField(-1, 'I')
1595 t
= z
.coefficient_tuple()
1600 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1601 [-c
+ d
*i
, a
- b
*i
]])
1602 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1603 blocks
.append(realM
)
1605 # We should have real entries by now, so use the realest field
1606 # we've got for the return value.
1607 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1612 def real_unembed(M
):
1614 The inverse of _embed_quaternion_matrix().
1618 sage: from mjo.eja.eja_algebra import \
1619 ....: QuaternionMatrixEuclideanJordanAlgebra
1623 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1624 ....: [-2, 1, -4, 3],
1625 ....: [-3, 4, 1, -2],
1626 ....: [-4, -3, 2, 1]])
1627 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1628 [1 + 2*i + 3*j + 4*k]
1632 Unembedding is the inverse of embedding::
1634 sage: set_random_seed()
1635 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1636 sage: M = random_matrix(Q, 3)
1637 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1638 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1644 raise ValueError("the matrix 'M' must be square")
1645 if not n
.mod(4).is_zero():
1646 raise ValueError("the matrix 'M' must be a quaternion embedding")
1648 # Use the base ring of the matrix to ensure that its entries can be
1649 # multiplied by elements of the quaternion algebra.
1650 field
= M
.base_ring()
1651 Q
= QuaternionAlgebra(field
,-1,-1)
1654 # Go top-left to bottom-right (reading order), converting every
1655 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1658 for l
in range(n
/4):
1659 for m
in range(n
/4):
1660 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1661 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1662 if submat
[0,0] != submat
[1,1].conjugate():
1663 raise ValueError('bad on-diagonal submatrix')
1664 if submat
[0,1] != -submat
[1,0].conjugate():
1665 raise ValueError('bad off-diagonal submatrix')
1666 z
= submat
[0,0].real()
1667 z
+= submat
[0,0].imag()*i
1668 z
+= submat
[0,1].real()*j
1669 z
+= submat
[0,1].imag()*k
1672 return matrix(Q
, n
/4, elements
)
1676 def natural_inner_product(cls
,X
,Y
):
1678 Compute a natural inner product in this algebra directly from
1683 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1687 This gives the same answer as the slow, default method implemented
1688 in :class:`MatrixEuclideanJordanAlgebra`::
1690 sage: set_random_seed()
1691 sage: J = QuaternionHermitianEJA.random_instance()
1692 sage: x,y = J.random_elements(2)
1693 sage: Xe = x.natural_representation()
1694 sage: Ye = y.natural_representation()
1695 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1696 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1697 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1698 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1699 sage: actual == expected
1703 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1706 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1708 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1709 matrices, the usual symmetric Jordan product, and the
1710 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1715 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1719 In theory, our "field" can be any subfield of the reals::
1721 sage: QuaternionHermitianEJA(2, RDF)
1722 Euclidean Jordan algebra of dimension 6 over Real Double Field
1723 sage: QuaternionHermitianEJA(2, RR)
1724 Euclidean Jordan algebra of dimension 6 over Real Field with
1725 53 bits of precision
1729 The dimension of this algebra is `2*n^2 - n`::
1731 sage: set_random_seed()
1732 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1733 sage: n = ZZ.random_element(1, n_max)
1734 sage: J = QuaternionHermitianEJA(n)
1735 sage: J.dimension() == 2*(n^2) - n
1738 The Jordan multiplication is what we think it is::
1740 sage: set_random_seed()
1741 sage: J = QuaternionHermitianEJA.random_instance()
1742 sage: x,y = J.random_elements(2)
1743 sage: actual = (x*y).natural_representation()
1744 sage: X = x.natural_representation()
1745 sage: Y = y.natural_representation()
1746 sage: expected = (X*Y + Y*X)/2
1747 sage: actual == expected
1749 sage: J(expected) == x*y
1752 We can change the generator prefix::
1754 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1755 (a0, a1, a2, a3, a4, a5)
1757 Our natural basis is normalized with respect to the natural inner
1758 product unless we specify otherwise::
1760 sage: set_random_seed()
1761 sage: J = QuaternionHermitianEJA.random_instance()
1762 sage: all( b.norm() == 1 for b in J.gens() )
1765 Since our natural basis is normalized with respect to the natural
1766 inner product, and since we know that this algebra is an EJA, any
1767 left-multiplication operator's matrix will be symmetric because
1768 natural->EJA basis representation is an isometry and within the EJA
1769 the operator is self-adjoint by the Jordan axiom::
1771 sage: set_random_seed()
1772 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1773 sage: x.operator().matrix().is_symmetric()
1776 We can construct the (trivial) algebra of rank zero::
1778 sage: QuaternionHermitianEJA(0)
1779 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1783 def _denormalized_basis(cls
, n
, field
):
1785 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1787 Why do we embed these? Basically, because all of numerical
1788 linear algebra assumes that you're working with vectors consisting
1789 of `n` entries from a field and scalars from the same field. There's
1790 no way to tell SageMath that (for example) the vectors contain
1791 complex numbers, while the scalar field is real.
1795 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1799 sage: set_random_seed()
1800 sage: n = ZZ.random_element(1,5)
1801 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1802 sage: all( M.is_symmetric() for M in B )
1806 Q
= QuaternionAlgebra(QQ
,-1,-1)
1809 # This is like the symmetric case, but we need to be careful:
1811 # * We want conjugate-symmetry, not just symmetry.
1812 # * The diagonal will (as a result) be real.
1816 for j
in range(i
+1):
1817 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1819 Sij
= cls
.real_embed(Eij
)
1822 # The second, third, and fourth ones have a minus
1823 # because they're conjugated.
1824 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1826 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1828 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1830 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1833 # Since we embedded these, we can drop back to the "field" that we
1834 # started with instead of the quaternion algebra "Q".
1835 return ( s
.change_ring(field
) for s
in S
)
1838 def __init__(self
, n
, field
=AA
, **kwargs
):
1839 basis
= self
._denormalized
_basis
(n
,field
)
1840 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1841 self
.rank
.set_cache(n
)
1844 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1846 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1847 with the half-trace inner product and jordan product ``x*y =
1848 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
1849 symmetric positive-definite "bilinear form" matrix. It has
1850 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
1851 when ``B`` is the identity matrix of order ``n-1``.
1855 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1856 ....: JordanSpinEJA)
1860 When no bilinear form is specified, the identity matrix is used,
1861 and the resulting algebra is the Jordan spin algebra::
1863 sage: J0 = BilinearFormEJA(3)
1864 sage: J1 = JordanSpinEJA(3)
1865 sage: J0.multiplication_table() == J0.multiplication_table()
1870 We can create a zero-dimensional algebra::
1872 sage: J = BilinearFormEJA(0)
1876 We can check the multiplication condition given in the Jordan, von
1877 Neumann, and Wigner paper (and also discussed on my "On the
1878 symmetry..." paper). Note that this relies heavily on the standard
1879 choice of basis, as does anything utilizing the bilinear form matrix::
1881 sage: set_random_seed()
1882 sage: n = ZZ.random_element(5)
1883 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
1884 sage: B = M.transpose()*M
1885 sage: J = BilinearFormEJA(n, B=B)
1886 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
1887 sage: V = J.vector_space()
1888 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
1889 ....: for ei in eis ]
1890 sage: actual = [ sis[i]*sis[j]
1891 ....: for i in range(n-1)
1892 ....: for j in range(n-1) ]
1893 sage: expected = [ J.one() if i == j else J.zero()
1894 ....: for i in range(n-1)
1895 ....: for j in range(n-1) ]
1896 sage: actual == expected
1899 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
1901 self
._B
= matrix
.identity(field
, max(0,n
-1))
1905 V
= VectorSpace(field
, n
)
1906 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1915 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
1916 zbar
= y0
*xbar
+ x0
*ybar
1917 z
= V([z0
] + zbar
.list())
1918 mult_table
[i
][j
] = z
1920 # The rank of this algebra is two, unless we're in a
1921 # one-dimensional ambient space (because the rank is bounded
1922 # by the ambient dimension).
1923 fdeja
= super(BilinearFormEJA
, self
)
1924 fdeja
.__init
__(field
, mult_table
, **kwargs
)
1925 self
.rank
.set_cache(min(n
,2))
1927 def inner_product(self
, x
, y
):
1929 Half of the trace inner product.
1931 This is defined so that the special case of the Jordan spin
1932 algebra gets the usual inner product.
1936 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1940 Ensure that this is one-half of the trace inner-product when
1941 the algebra isn't just the reals (when ``n`` isn't one). This
1942 is in Faraut and Koranyi, and also my "On the symmetry..."
1945 sage: set_random_seed()
1946 sage: n = ZZ.random_element(2,5)
1947 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
1948 sage: B = M.transpose()*M
1949 sage: J = BilinearFormEJA(n, B=B)
1950 sage: x = J.random_element()
1951 sage: y = J.random_element()
1952 sage: x.inner_product(y) == (x*y).trace()/2
1956 xvec
= x
.to_vector()
1958 yvec
= y
.to_vector()
1960 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
1963 class JordanSpinEJA(BilinearFormEJA
):
1965 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1966 with the usual inner product and jordan product ``x*y =
1967 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1972 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1976 This multiplication table can be verified by hand::
1978 sage: J = JordanSpinEJA(4)
1979 sage: e0,e1,e2,e3 = J.gens()
1995 We can change the generator prefix::
1997 sage: JordanSpinEJA(2, prefix='B').gens()
2002 Ensure that we have the usual inner product on `R^n`::
2004 sage: set_random_seed()
2005 sage: J = JordanSpinEJA.random_instance()
2006 sage: x,y = J.random_elements(2)
2007 sage: X = x.natural_representation()
2008 sage: Y = y.natural_representation()
2009 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2013 def __init__(self
, n
, field
=AA
, **kwargs
):
2014 # This is a special case of the BilinearFormEJA with the identity
2015 # matrix as its bilinear form.
2016 return super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2019 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2021 The trivial Euclidean Jordan algebra consisting of only a zero element.
2025 sage: from mjo.eja.eja_algebra import TrivialEJA
2029 sage: J = TrivialEJA()
2036 sage: 7*J.one()*12*J.one()
2038 sage: J.one().inner_product(J.one())
2040 sage: J.one().norm()
2042 sage: J.one().subalgebra_generated_by()
2043 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2048 def __init__(self
, field
=AA
, **kwargs
):
2050 fdeja
= super(TrivialEJA
, self
)
2051 # The rank is zero using my definition, namely the dimension of the
2052 # largest subalgebra generated by any element.
2053 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2054 self
.rank
.set_cache(0)