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...
584 sage: J0.one().natural_representation()
588 sage: orig_df = AA.options.display_format
589 sage: AA.options.display_format = 'radical'
590 sage: J.from_vector(J5.basis()[0]).natural_representation()
594 sage: J.from_vector(J5.basis()[1]).natural_representation()
598 sage: AA.options.display_format = orig_df
599 sage: J1.one().natural_representation()
606 Every algebra decomposes trivially with respect to its identity
609 sage: set_random_seed()
610 sage: J = random_eja()
611 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
612 sage: J0.dimension() == 0 and J5.dimension() == 0
614 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
617 The decomposition is into eigenspaces, and its components are
618 therefore necessarily orthogonal. Moreover, the identity
619 elements in the two subalgebras are the projections onto their
620 respective subspaces of the superalgebra's identity element::
622 sage: set_random_seed()
623 sage: J = random_eja()
624 sage: x = J.random_element()
625 sage: if not J.is_trivial():
626 ....: while x.is_nilpotent():
627 ....: x = J.random_element()
628 sage: c = x.subalgebra_idempotent()
629 sage: J0,J5,J1 = J.peirce_decomposition(c)
631 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
632 ....: w = w.superalgebra_element()
633 ....: y = J.from_vector(y)
634 ....: z = z.superalgebra_element()
635 ....: ipsum += w.inner_product(y).abs()
636 ....: ipsum += w.inner_product(z).abs()
637 ....: ipsum += y.inner_product(z).abs()
640 sage: J1(c) == J1.one()
642 sage: J0(J.one() - c) == J0.one()
646 if not c
.is_idempotent():
647 raise ValueError("element is not idempotent: %s" % c
)
649 # Default these to what they should be if they turn out to be
650 # trivial, because eigenspaces_left() won't return eigenvalues
651 # corresponding to trivial spaces (e.g. it returns only the
652 # eigenspace corresponding to lambda=1 if you take the
653 # decomposition relative to the identity element).
654 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
655 J0
= trivial
# eigenvalue zero
656 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
657 J1
= trivial
# eigenvalue one
659 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
660 if eigval
== ~
(self
.base_ring()(2)):
663 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
664 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
670 raise ValueError("unexpected eigenvalue: %s" % eigval
)
675 def random_element(self
, thorough
=False):
677 Return a random element of this algebra.
679 Our algebra superclass method only returns a linear
680 combination of at most two basis elements. We instead
681 want the vector space "random element" method that
682 returns a more diverse selection.
686 - ``thorough`` -- (boolean; default False) whether or not we
687 should generate irrational coefficients for the random
688 element when our base ring is irrational; this slows the
689 algebra operations to a crawl, but any truly random method
693 # For a general base ring... maybe we can trust this to do the
694 # right thing? Unlikely, but.
695 V
= self
.vector_space()
696 v
= V
.random_element()
698 if self
.base_ring() is AA
:
699 # The "random element" method of the algebraic reals is
700 # stupid at the moment, and only returns integers between
701 # -2 and 2, inclusive:
703 # https://trac.sagemath.org/ticket/30875
705 # Instead, we implement our own "random vector" method,
706 # and then coerce that into the algebra. We use the vector
707 # space degree here instead of the dimension because a
708 # subalgebra could (for example) be spanned by only two
709 # vectors, each with five coordinates. We need to
710 # generate all five coordinates.
712 v
*= QQbar
.random_element().real()
714 v
*= QQ
.random_element()
716 return self
.from_vector(V
.coordinate_vector(v
))
718 def random_elements(self
, count
, thorough
=False):
720 Return ``count`` random elements as a tuple.
724 - ``thorough`` -- (boolean; default False) whether or not we
725 should generate irrational coefficients for the random
726 elements when our base ring is irrational; this slows the
727 algebra operations to a crawl, but any truly random method
732 sage: from mjo.eja.eja_algebra import JordanSpinEJA
736 sage: J = JordanSpinEJA(3)
737 sage: x,y,z = J.random_elements(3)
738 sage: all( [ x in J, y in J, z in J ])
740 sage: len( J.random_elements(10) ) == 10
744 return tuple( self
.random_element(thorough
)
745 for idx
in range(count
) )
748 def random_instance(cls
, field
=AA
, **kwargs
):
750 Return a random instance of this type of algebra.
752 Beware, this will crash for "most instances" because the
753 constructor below looks wrong.
755 if cls
is TrivialEJA
:
756 # The TrivialEJA class doesn't take an "n" argument because
760 n
= ZZ
.random_element(cls
._max
_test
_case
_size
() + 1)
761 return cls(n
, field
, **kwargs
)
764 def _charpoly_coefficients(self
):
766 The `r` polynomial coefficients of the "characteristic polynomial
770 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
771 R
= PolynomialRing(self
.base_ring(), var_names
)
773 F
= R
.fraction_field()
776 # From a result in my book, these are the entries of the
777 # basis representation of L_x.
778 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
781 L_x
= matrix(F
, n
, n
, L_x_i_j
)
784 if self
.rank
.is_in_cache():
786 # There's no need to pad the system with redundant
787 # columns if we *know* they'll be redundant.
790 # Compute an extra power in case the rank is equal to
791 # the dimension (otherwise, we would stop at x^(r-1)).
792 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
793 for k
in range(n
+1) ]
794 A
= matrix
.column(F
, x_powers
[:n
])
795 AE
= A
.extended_echelon_form()
802 # The theory says that only the first "r" coefficients are
803 # nonzero, and they actually live in the original polynomial
804 # ring and not the fraction field. We negate them because
805 # in the actual characteristic polynomial, they get moved
806 # to the other side where x^r lives.
807 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
812 Return the rank of this EJA.
814 This is a cached method because we know the rank a priori for
815 all of the algebras we can construct. Thus we can avoid the
816 expensive ``_charpoly_coefficients()`` call unless we truly
817 need to compute the whole characteristic polynomial.
821 sage: from mjo.eja.eja_algebra import (HadamardEJA,
823 ....: RealSymmetricEJA,
824 ....: ComplexHermitianEJA,
825 ....: QuaternionHermitianEJA,
830 The rank of the Jordan spin algebra is always two::
832 sage: JordanSpinEJA(2).rank()
834 sage: JordanSpinEJA(3).rank()
836 sage: JordanSpinEJA(4).rank()
839 The rank of the `n`-by-`n` Hermitian real, complex, or
840 quaternion matrices is `n`::
842 sage: RealSymmetricEJA(4).rank()
844 sage: ComplexHermitianEJA(3).rank()
846 sage: QuaternionHermitianEJA(2).rank()
851 Ensure that every EJA that we know how to construct has a
852 positive integer rank, unless the algebra is trivial in
853 which case its rank will be zero::
855 sage: set_random_seed()
856 sage: J = random_eja()
860 sage: r > 0 or (r == 0 and J.is_trivial())
863 Ensure that computing the rank actually works, since the ranks
864 of all simple algebras are known and will be cached by default::
866 sage: J = HadamardEJA(4)
867 sage: J.rank.clear_cache()
873 sage: J = JordanSpinEJA(4)
874 sage: J.rank.clear_cache()
880 sage: J = RealSymmetricEJA(3)
881 sage: J.rank.clear_cache()
887 sage: J = ComplexHermitianEJA(2)
888 sage: J.rank.clear_cache()
894 sage: J = QuaternionHermitianEJA(2)
895 sage: J.rank.clear_cache()
899 return len(self
._charpoly
_coefficients
())
902 def vector_space(self
):
904 Return the vector space that underlies this algebra.
908 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
912 sage: J = RealSymmetricEJA(2)
913 sage: J.vector_space()
914 Vector space of dimension 3 over...
917 return self
.zero().to_vector().parent().ambient_vector_space()
920 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
923 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
):
925 Return the Euclidean Jordan Algebra corresponding to the set
926 `R^n` under the Hadamard product.
928 Note: this is nothing more than the Cartesian product of ``n``
929 copies of the spin algebra. Once Cartesian product algebras
930 are implemented, this can go.
934 sage: from mjo.eja.eja_algebra import HadamardEJA
938 This multiplication table can be verified by hand::
940 sage: J = HadamardEJA(3)
941 sage: e0,e1,e2 = J.gens()
957 We can change the generator prefix::
959 sage: HadamardEJA(3, prefix='r').gens()
963 def __init__(self
, n
, field
=AA
, **kwargs
):
964 V
= VectorSpace(field
, n
)
965 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
968 fdeja
= super(HadamardEJA
, self
)
969 fdeja
.__init
__(field
, mult_table
, **kwargs
)
970 self
.rank
.set_cache(n
)
972 def inner_product(self
, x
, y
):
974 Faster to reimplement than to use natural representations.
978 sage: from mjo.eja.eja_algebra import HadamardEJA
982 Ensure that this is the usual inner product for the algebras
985 sage: set_random_seed()
986 sage: J = HadamardEJA.random_instance()
987 sage: x,y = J.random_elements(2)
988 sage: X = x.natural_representation()
989 sage: Y = y.natural_representation()
990 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
994 return x
.to_vector().inner_product(y
.to_vector())
997 def random_eja(field
=AA
):
999 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1003 sage: from mjo.eja.eja_algebra import random_eja
1008 Euclidean Jordan algebra of dimension...
1011 classname
= choice([TrivialEJA
,
1015 ComplexHermitianEJA
,
1016 QuaternionHermitianEJA
])
1017 return classname
.random_instance(field
=field
)
1022 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1024 def _max_test_case_size():
1025 # Play it safe, since this will be squared and the underlying
1026 # field can have dimension 4 (quaternions) too.
1029 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1031 Compared to the superclass constructor, we take a basis instead of
1032 a multiplication table because the latter can be computed in terms
1033 of the former when the product is known (like it is here).
1035 # Used in this class's fast _charpoly_coefficients() override.
1036 self
._basis
_normalizers
= None
1038 # We're going to loop through this a few times, so now's a good
1039 # time to ensure that it isn't a generator expression.
1040 basis
= tuple(basis
)
1042 if len(basis
) > 1 and normalize_basis
:
1043 # We'll need sqrt(2) to normalize the basis, and this
1044 # winds up in the multiplication table, so the whole
1045 # algebra needs to be over the field extension.
1046 R
= PolynomialRing(field
, 'z')
1049 if p
.is_irreducible():
1050 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1051 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1052 self
._basis
_normalizers
= tuple(
1053 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1054 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1056 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1058 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
1059 fdeja
.__init
__(field
, Qs
, natural_basis
=basis
, **kwargs
)
1064 def _charpoly_coefficients(self
):
1066 Override the parent method with something that tries to compute
1067 over a faster (non-extension) field.
1069 if self
._basis
_normalizers
is None:
1070 # We didn't normalize, so assume that the basis we started
1071 # with had entries in a nice field.
1072 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1074 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1075 self
._basis
_normalizers
) )
1077 # Do this over the rationals and convert back at the end.
1078 # Only works because we know the entries of the basis are
1080 J
= MatrixEuclideanJordanAlgebra(QQ
,
1082 normalize_basis
=False)
1083 a
= J
._charpoly
_coefficients
()
1085 # Unfortunately, changing the basis does change the
1086 # coefficients of the characteristic polynomial, but since
1087 # these are really the coefficients of the "characteristic
1088 # polynomial of" function, everything is still nice and
1089 # unevaluated. It's therefore "obvious" how scaling the
1090 # basis affects the coordinate variables X1, X2, et
1091 # cetera. Scaling the first basis vector up by "n" adds a
1092 # factor of 1/n into every "X1" term, for example. So here
1093 # we simply undo the basis_normalizer scaling that we
1094 # performed earlier.
1096 # The a[0] access here is safe because trivial algebras
1097 # won't have any basis normalizers and therefore won't
1098 # make it to this "else" branch.
1099 XS
= a
[0].parent().gens()
1100 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1101 for i
in range(len(XS
)) }
1102 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1106 def multiplication_table_from_matrix_basis(basis
):
1108 At least three of the five simple Euclidean Jordan algebras have the
1109 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1110 multiplication on the right is matrix multiplication. Given a basis
1111 for the underlying matrix space, this function returns a
1112 multiplication table (obtained by looping through the basis
1113 elements) for an algebra of those matrices.
1115 # In S^2, for example, we nominally have four coordinates even
1116 # though the space is of dimension three only. The vector space V
1117 # is supposed to hold the entire long vector, and the subspace W
1118 # of V will be spanned by the vectors that arise from symmetric
1119 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1123 field
= basis
[0].base_ring()
1124 dimension
= basis
[0].nrows()
1126 V
= VectorSpace(field
, dimension
**2)
1127 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1129 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1132 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1133 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1141 Embed the matrix ``M`` into a space of real matrices.
1143 The matrix ``M`` can have entries in any field at the moment:
1144 the real numbers, complex numbers, or quaternions. And although
1145 they are not a field, we can probably support octonions at some
1146 point, too. This function returns a real matrix that "acts like"
1147 the original with respect to matrix multiplication; i.e.
1149 real_embed(M*N) = real_embed(M)*real_embed(N)
1152 raise NotImplementedError
1156 def real_unembed(M
):
1158 The inverse of :meth:`real_embed`.
1160 raise NotImplementedError
1164 def natural_inner_product(cls
,X
,Y
):
1165 Xu
= cls
.real_unembed(X
)
1166 Yu
= cls
.real_unembed(Y
)
1167 tr
= (Xu
*Yu
).trace()
1170 # It's real already.
1173 # Otherwise, try the thing that works for complex numbers; and
1174 # if that doesn't work, the thing that works for quaternions.
1176 return tr
.vector()[0] # real part, imag part is index 1
1177 except AttributeError:
1178 # A quaternions doesn't have a vector() method, but does
1179 # have coefficient_tuple() method that returns the
1180 # coefficients of 1, i, j, and k -- in that order.
1181 return tr
.coefficient_tuple()[0]
1184 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1188 The identity function, for embedding real matrices into real
1194 def real_unembed(M
):
1196 The identity function, for unembedding real matrices from real
1202 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1204 The rank-n simple EJA consisting of real symmetric n-by-n
1205 matrices, the usual symmetric Jordan product, and the trace inner
1206 product. It has dimension `(n^2 + n)/2` over the reals.
1210 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1214 sage: J = RealSymmetricEJA(2)
1215 sage: e0, e1, e2 = J.gens()
1223 In theory, our "field" can be any subfield of the reals::
1225 sage: RealSymmetricEJA(2, RDF)
1226 Euclidean Jordan algebra of dimension 3 over Real Double Field
1227 sage: RealSymmetricEJA(2, RR)
1228 Euclidean Jordan algebra of dimension 3 over Real Field with
1229 53 bits of precision
1233 The dimension of this algebra is `(n^2 + n) / 2`::
1235 sage: set_random_seed()
1236 sage: n_max = RealSymmetricEJA._max_test_case_size()
1237 sage: n = ZZ.random_element(1, n_max)
1238 sage: J = RealSymmetricEJA(n)
1239 sage: J.dimension() == (n^2 + n)/2
1242 The Jordan multiplication is what we think it is::
1244 sage: set_random_seed()
1245 sage: J = RealSymmetricEJA.random_instance()
1246 sage: x,y = J.random_elements(2)
1247 sage: actual = (x*y).natural_representation()
1248 sage: X = x.natural_representation()
1249 sage: Y = y.natural_representation()
1250 sage: expected = (X*Y + Y*X)/2
1251 sage: actual == expected
1253 sage: J(expected) == x*y
1256 We can change the generator prefix::
1258 sage: RealSymmetricEJA(3, prefix='q').gens()
1259 (q0, q1, q2, q3, q4, q5)
1261 Our natural basis is normalized with respect to the natural inner
1262 product unless we specify otherwise::
1264 sage: set_random_seed()
1265 sage: J = RealSymmetricEJA.random_instance()
1266 sage: all( b.norm() == 1 for b in J.gens() )
1269 Since our natural basis is normalized with respect to the natural
1270 inner product, and since we know that this algebra is an EJA, any
1271 left-multiplication operator's matrix will be symmetric because
1272 natural->EJA basis representation is an isometry and within the EJA
1273 the operator is self-adjoint by the Jordan axiom::
1275 sage: set_random_seed()
1276 sage: x = RealSymmetricEJA.random_instance().random_element()
1277 sage: x.operator().matrix().is_symmetric()
1280 We can construct the (trivial) algebra of rank zero::
1282 sage: RealSymmetricEJA(0)
1283 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1287 def _denormalized_basis(cls
, n
, field
):
1289 Return a basis for the space of real symmetric n-by-n matrices.
1293 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1297 sage: set_random_seed()
1298 sage: n = ZZ.random_element(1,5)
1299 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1300 sage: all( M.is_symmetric() for M in B)
1304 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1308 for j
in range(i
+1):
1309 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1313 Sij
= Eij
+ Eij
.transpose()
1319 def _max_test_case_size():
1320 return 4 # Dimension 10
1323 def __init__(self
, n
, field
=AA
, **kwargs
):
1324 basis
= self
._denormalized
_basis
(n
, field
)
1325 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, **kwargs
)
1326 self
.rank
.set_cache(n
)
1329 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1333 Embed the n-by-n complex matrix ``M`` into the space of real
1334 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1335 bi` to the block matrix ``[[a,b],[-b,a]]``.
1339 sage: from mjo.eja.eja_algebra import \
1340 ....: ComplexMatrixEuclideanJordanAlgebra
1344 sage: F = QuadraticField(-1, 'I')
1345 sage: x1 = F(4 - 2*i)
1346 sage: x2 = F(1 + 2*i)
1349 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1350 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1359 Embedding is a homomorphism (isomorphism, in fact)::
1361 sage: set_random_seed()
1362 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1363 sage: n = ZZ.random_element(n_max)
1364 sage: F = QuadraticField(-1, 'I')
1365 sage: X = random_matrix(F, n)
1366 sage: Y = random_matrix(F, n)
1367 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1368 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1369 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1376 raise ValueError("the matrix 'M' must be square")
1378 # We don't need any adjoined elements...
1379 field
= M
.base_ring().base_ring()
1383 a
= z
.list()[0] # real part, I guess
1384 b
= z
.list()[1] # imag part, I guess
1385 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1387 return matrix
.block(field
, n
, blocks
)
1391 def real_unembed(M
):
1393 The inverse of _embed_complex_matrix().
1397 sage: from mjo.eja.eja_algebra import \
1398 ....: ComplexMatrixEuclideanJordanAlgebra
1402 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1403 ....: [-2, 1, -4, 3],
1404 ....: [ 9, 10, 11, 12],
1405 ....: [-10, 9, -12, 11] ])
1406 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1408 [ 10*I + 9 12*I + 11]
1412 Unembedding is the inverse of embedding::
1414 sage: set_random_seed()
1415 sage: F = QuadraticField(-1, 'I')
1416 sage: M = random_matrix(F, 3)
1417 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1418 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1424 raise ValueError("the matrix 'M' must be square")
1425 if not n
.mod(2).is_zero():
1426 raise ValueError("the matrix 'M' must be a complex embedding")
1428 # If "M" was normalized, its base ring might have roots
1429 # adjoined and they can stick around after unembedding.
1430 field
= M
.base_ring()
1431 R
= PolynomialRing(field
, 'z')
1434 # Sage doesn't know how to embed AA into QQbar, i.e. how
1435 # to adjoin sqrt(-1) to AA.
1438 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1441 # Go top-left to bottom-right (reading order), converting every
1442 # 2-by-2 block we see to a single complex element.
1444 for k
in range(n
/2):
1445 for j
in range(n
/2):
1446 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1447 if submat
[0,0] != submat
[1,1]:
1448 raise ValueError('bad on-diagonal submatrix')
1449 if submat
[0,1] != -submat
[1,0]:
1450 raise ValueError('bad off-diagonal submatrix')
1451 z
= submat
[0,0] + submat
[0,1]*i
1454 return matrix(F
, n
/2, elements
)
1458 def natural_inner_product(cls
,X
,Y
):
1460 Compute a natural inner product in this algebra directly from
1465 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1469 This gives the same answer as the slow, default method implemented
1470 in :class:`MatrixEuclideanJordanAlgebra`::
1472 sage: set_random_seed()
1473 sage: J = ComplexHermitianEJA.random_instance()
1474 sage: x,y = J.random_elements(2)
1475 sage: Xe = x.natural_representation()
1476 sage: Ye = y.natural_representation()
1477 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1478 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1479 sage: expected = (X*Y).trace().real()
1480 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1481 sage: actual == expected
1485 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1488 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1490 The rank-n simple EJA consisting of complex Hermitian n-by-n
1491 matrices over the real numbers, the usual symmetric Jordan product,
1492 and the real-part-of-trace inner product. It has dimension `n^2` over
1497 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1501 In theory, our "field" can be any subfield of the reals::
1503 sage: ComplexHermitianEJA(2, RDF)
1504 Euclidean Jordan algebra of dimension 4 over Real Double Field
1505 sage: ComplexHermitianEJA(2, RR)
1506 Euclidean Jordan algebra of dimension 4 over Real Field with
1507 53 bits of precision
1511 The dimension of this algebra is `n^2`::
1513 sage: set_random_seed()
1514 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1515 sage: n = ZZ.random_element(1, n_max)
1516 sage: J = ComplexHermitianEJA(n)
1517 sage: J.dimension() == n^2
1520 The Jordan multiplication is what we think it is::
1522 sage: set_random_seed()
1523 sage: J = ComplexHermitianEJA.random_instance()
1524 sage: x,y = J.random_elements(2)
1525 sage: actual = (x*y).natural_representation()
1526 sage: X = x.natural_representation()
1527 sage: Y = y.natural_representation()
1528 sage: expected = (X*Y + Y*X)/2
1529 sage: actual == expected
1531 sage: J(expected) == x*y
1534 We can change the generator prefix::
1536 sage: ComplexHermitianEJA(2, prefix='z').gens()
1539 Our natural basis is normalized with respect to the natural inner
1540 product unless we specify otherwise::
1542 sage: set_random_seed()
1543 sage: J = ComplexHermitianEJA.random_instance()
1544 sage: all( b.norm() == 1 for b in J.gens() )
1547 Since our natural basis is normalized with respect to the natural
1548 inner product, and since we know that this algebra is an EJA, any
1549 left-multiplication operator's matrix will be symmetric because
1550 natural->EJA basis representation is an isometry and within the EJA
1551 the operator is self-adjoint by the Jordan axiom::
1553 sage: set_random_seed()
1554 sage: x = ComplexHermitianEJA.random_instance().random_element()
1555 sage: x.operator().matrix().is_symmetric()
1558 We can construct the (trivial) algebra of rank zero::
1560 sage: ComplexHermitianEJA(0)
1561 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1566 def _denormalized_basis(cls
, n
, field
):
1568 Returns a basis for the space of complex Hermitian n-by-n matrices.
1570 Why do we embed these? Basically, because all of numerical linear
1571 algebra assumes that you're working with vectors consisting of `n`
1572 entries from a field and scalars from the same field. There's no way
1573 to tell SageMath that (for example) the vectors contain complex
1574 numbers, while the scalar field is real.
1578 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1582 sage: set_random_seed()
1583 sage: n = ZZ.random_element(1,5)
1584 sage: field = QuadraticField(2, 'sqrt2')
1585 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1586 sage: all( M.is_symmetric() for M in B)
1590 R
= PolynomialRing(field
, 'z')
1592 F
= field
.extension(z
**2 + 1, 'I')
1595 # This is like the symmetric case, but we need to be careful:
1597 # * We want conjugate-symmetry, not just symmetry.
1598 # * The diagonal will (as a result) be real.
1602 for j
in range(i
+1):
1603 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1605 Sij
= cls
.real_embed(Eij
)
1608 # The second one has a minus because it's conjugated.
1609 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1611 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1614 # Since we embedded these, we can drop back to the "field" that we
1615 # started with instead of the complex extension "F".
1616 return ( s
.change_ring(field
) for s
in S
)
1619 def __init__(self
, n
, field
=AA
, **kwargs
):
1620 basis
= self
._denormalized
_basis
(n
,field
)
1621 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1622 self
.rank
.set_cache(n
)
1625 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1629 Embed the n-by-n quaternion matrix ``M`` into the space of real
1630 matrices of size 4n-by-4n by first sending each quaternion entry `z
1631 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1632 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1637 sage: from mjo.eja.eja_algebra import \
1638 ....: QuaternionMatrixEuclideanJordanAlgebra
1642 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1643 sage: i,j,k = Q.gens()
1644 sage: x = 1 + 2*i + 3*j + 4*k
1645 sage: M = matrix(Q, 1, [[x]])
1646 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1652 Embedding is a homomorphism (isomorphism, in fact)::
1654 sage: set_random_seed()
1655 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1656 sage: n = ZZ.random_element(n_max)
1657 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1658 sage: X = random_matrix(Q, n)
1659 sage: Y = random_matrix(Q, n)
1660 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1661 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1662 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1667 quaternions
= M
.base_ring()
1670 raise ValueError("the matrix 'M' must be square")
1672 F
= QuadraticField(-1, 'I')
1677 t
= z
.coefficient_tuple()
1682 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1683 [-c
+ d
*i
, a
- b
*i
]])
1684 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1685 blocks
.append(realM
)
1687 # We should have real entries by now, so use the realest field
1688 # we've got for the return value.
1689 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1694 def real_unembed(M
):
1696 The inverse of _embed_quaternion_matrix().
1700 sage: from mjo.eja.eja_algebra import \
1701 ....: QuaternionMatrixEuclideanJordanAlgebra
1705 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1706 ....: [-2, 1, -4, 3],
1707 ....: [-3, 4, 1, -2],
1708 ....: [-4, -3, 2, 1]])
1709 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1710 [1 + 2*i + 3*j + 4*k]
1714 Unembedding is the inverse of embedding::
1716 sage: set_random_seed()
1717 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1718 sage: M = random_matrix(Q, 3)
1719 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1720 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1726 raise ValueError("the matrix 'M' must be square")
1727 if not n
.mod(4).is_zero():
1728 raise ValueError("the matrix 'M' must be a quaternion embedding")
1730 # Use the base ring of the matrix to ensure that its entries can be
1731 # multiplied by elements of the quaternion algebra.
1732 field
= M
.base_ring()
1733 Q
= QuaternionAlgebra(field
,-1,-1)
1736 # Go top-left to bottom-right (reading order), converting every
1737 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1740 for l
in range(n
/4):
1741 for m
in range(n
/4):
1742 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1743 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1744 if submat
[0,0] != submat
[1,1].conjugate():
1745 raise ValueError('bad on-diagonal submatrix')
1746 if submat
[0,1] != -submat
[1,0].conjugate():
1747 raise ValueError('bad off-diagonal submatrix')
1748 z
= submat
[0,0].real()
1749 z
+= submat
[0,0].imag()*i
1750 z
+= submat
[0,1].real()*j
1751 z
+= submat
[0,1].imag()*k
1754 return matrix(Q
, n
/4, elements
)
1758 def natural_inner_product(cls
,X
,Y
):
1760 Compute a natural inner product in this algebra directly from
1765 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1769 This gives the same answer as the slow, default method implemented
1770 in :class:`MatrixEuclideanJordanAlgebra`::
1772 sage: set_random_seed()
1773 sage: J = QuaternionHermitianEJA.random_instance()
1774 sage: x,y = J.random_elements(2)
1775 sage: Xe = x.natural_representation()
1776 sage: Ye = y.natural_representation()
1777 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1778 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1779 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1780 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1781 sage: actual == expected
1785 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1788 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1790 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1791 matrices, the usual symmetric Jordan product, and the
1792 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1797 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1801 In theory, our "field" can be any subfield of the reals::
1803 sage: QuaternionHermitianEJA(2, RDF)
1804 Euclidean Jordan algebra of dimension 6 over Real Double Field
1805 sage: QuaternionHermitianEJA(2, RR)
1806 Euclidean Jordan algebra of dimension 6 over Real Field with
1807 53 bits of precision
1811 The dimension of this algebra is `2*n^2 - n`::
1813 sage: set_random_seed()
1814 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1815 sage: n = ZZ.random_element(1, n_max)
1816 sage: J = QuaternionHermitianEJA(n)
1817 sage: J.dimension() == 2*(n^2) - n
1820 The Jordan multiplication is what we think it is::
1822 sage: set_random_seed()
1823 sage: J = QuaternionHermitianEJA.random_instance()
1824 sage: x,y = J.random_elements(2)
1825 sage: actual = (x*y).natural_representation()
1826 sage: X = x.natural_representation()
1827 sage: Y = y.natural_representation()
1828 sage: expected = (X*Y + Y*X)/2
1829 sage: actual == expected
1831 sage: J(expected) == x*y
1834 We can change the generator prefix::
1836 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1837 (a0, a1, a2, a3, a4, a5)
1839 Our natural basis is normalized with respect to the natural inner
1840 product unless we specify otherwise::
1842 sage: set_random_seed()
1843 sage: J = QuaternionHermitianEJA.random_instance()
1844 sage: all( b.norm() == 1 for b in J.gens() )
1847 Since our natural basis is normalized with respect to the natural
1848 inner product, and since we know that this algebra is an EJA, any
1849 left-multiplication operator's matrix will be symmetric because
1850 natural->EJA basis representation is an isometry and within the EJA
1851 the operator is self-adjoint by the Jordan axiom::
1853 sage: set_random_seed()
1854 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1855 sage: x.operator().matrix().is_symmetric()
1858 We can construct the (trivial) algebra of rank zero::
1860 sage: QuaternionHermitianEJA(0)
1861 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1865 def _denormalized_basis(cls
, n
, field
):
1867 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1869 Why do we embed these? Basically, because all of numerical
1870 linear algebra assumes that you're working with vectors consisting
1871 of `n` entries from a field and scalars from the same field. There's
1872 no way to tell SageMath that (for example) the vectors contain
1873 complex numbers, while the scalar field is real.
1877 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1881 sage: set_random_seed()
1882 sage: n = ZZ.random_element(1,5)
1883 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1884 sage: all( M.is_symmetric() for M in B )
1888 Q
= QuaternionAlgebra(QQ
,-1,-1)
1891 # This is like the symmetric case, but we need to be careful:
1893 # * We want conjugate-symmetry, not just symmetry.
1894 # * The diagonal will (as a result) be real.
1898 for j
in range(i
+1):
1899 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1901 Sij
= cls
.real_embed(Eij
)
1904 # The second, third, and fourth ones have a minus
1905 # because they're conjugated.
1906 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1908 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1910 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1912 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1915 # Since we embedded these, we can drop back to the "field" that we
1916 # started with instead of the quaternion algebra "Q".
1917 return ( s
.change_ring(field
) for s
in S
)
1920 def __init__(self
, n
, field
=AA
, **kwargs
):
1921 basis
= self
._denormalized
_basis
(n
,field
)
1922 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1923 self
.rank
.set_cache(n
)
1926 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1928 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1929 with the half-trace inner product and jordan product ``x*y =
1930 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
1931 symmetric positive-definite "bilinear form" matrix. It has
1932 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
1933 when ``B`` is the identity matrix of order ``n-1``.
1937 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1938 ....: JordanSpinEJA)
1942 When no bilinear form is specified, the identity matrix is used,
1943 and the resulting algebra is the Jordan spin algebra::
1945 sage: J0 = BilinearFormEJA(3)
1946 sage: J1 = JordanSpinEJA(3)
1947 sage: J0.multiplication_table() == J0.multiplication_table()
1952 We can create a zero-dimensional algebra::
1954 sage: J = BilinearFormEJA(0)
1958 We can check the multiplication condition given in the Jordan, von
1959 Neumann, and Wigner paper (and also discussed on my "On the
1960 symmetry..." paper). Note that this relies heavily on the standard
1961 choice of basis, as does anything utilizing the bilinear form matrix::
1963 sage: set_random_seed()
1964 sage: n = ZZ.random_element(5)
1965 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
1966 sage: B = M.transpose()*M
1967 sage: J = BilinearFormEJA(n, B=B)
1968 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
1969 sage: V = J.vector_space()
1970 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
1971 ....: for ei in eis ]
1972 sage: actual = [ sis[i]*sis[j]
1973 ....: for i in range(n-1)
1974 ....: for j in range(n-1) ]
1975 sage: expected = [ J.one() if i == j else J.zero()
1976 ....: for i in range(n-1)
1977 ....: for j in range(n-1) ]
1978 sage: actual == expected
1981 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
1983 self
._B
= matrix
.identity(field
, max(0,n
-1))
1987 V
= VectorSpace(field
, n
)
1988 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1997 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
1998 zbar
= y0
*xbar
+ x0
*ybar
1999 z
= V([z0
] + zbar
.list())
2000 mult_table
[i
][j
] = z
2002 # The rank of this algebra is two, unless we're in a
2003 # one-dimensional ambient space (because the rank is bounded
2004 # by the ambient dimension).
2005 fdeja
= super(BilinearFormEJA
, self
)
2006 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2007 self
.rank
.set_cache(min(n
,2))
2009 def inner_product(self
, x
, y
):
2011 Half of the trace inner product.
2013 This is defined so that the special case of the Jordan spin
2014 algebra gets the usual inner product.
2018 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2022 Ensure that this is one-half of the trace inner-product when
2023 the algebra isn't just the reals (when ``n`` isn't one). This
2024 is in Faraut and Koranyi, and also my "On the symmetry..."
2027 sage: set_random_seed()
2028 sage: n = ZZ.random_element(2,5)
2029 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2030 sage: B = M.transpose()*M
2031 sage: J = BilinearFormEJA(n, B=B)
2032 sage: x = J.random_element()
2033 sage: y = J.random_element()
2034 sage: x.inner_product(y) == (x*y).trace()/2
2038 xvec
= x
.to_vector()
2040 yvec
= y
.to_vector()
2042 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
2045 class JordanSpinEJA(BilinearFormEJA
):
2047 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2048 with the usual inner product and jordan product ``x*y =
2049 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2054 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2058 This multiplication table can be verified by hand::
2060 sage: J = JordanSpinEJA(4)
2061 sage: e0,e1,e2,e3 = J.gens()
2077 We can change the generator prefix::
2079 sage: JordanSpinEJA(2, prefix='B').gens()
2084 Ensure that we have the usual inner product on `R^n`::
2086 sage: set_random_seed()
2087 sage: J = JordanSpinEJA.random_instance()
2088 sage: x,y = J.random_elements(2)
2089 sage: X = x.natural_representation()
2090 sage: Y = y.natural_representation()
2091 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2095 def __init__(self
, n
, field
=AA
, **kwargs
):
2096 # This is a special case of the BilinearFormEJA with the identity
2097 # matrix as its bilinear form.
2098 return super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2101 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2103 The trivial Euclidean Jordan algebra consisting of only a zero element.
2107 sage: from mjo.eja.eja_algebra import TrivialEJA
2111 sage: J = TrivialEJA()
2118 sage: 7*J.one()*12*J.one()
2120 sage: J.one().inner_product(J.one())
2122 sage: J.one().norm()
2124 sage: J.one().subalgebra_generated_by()
2125 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2130 def __init__(self
, field
=AA
, **kwargs
):
2132 fdeja
= super(TrivialEJA
, self
)
2133 # The rank is zero using my definition, namely the dimension of the
2134 # largest subalgebra generated by any element.
2135 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2136 self
.rank
.set_cache(0)
2139 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2141 The external (orthogonal) direct sum of two other Euclidean Jordan
2142 algebras. Essentially the Cartesian product of its two factors.
2143 Every Euclidean Jordan algebra decomposes into an orthogonal
2144 direct sum of simple Euclidean Jordan algebras, so no generality
2145 is lost by providing only this construction.
2149 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2150 ....: RealSymmetricEJA,
2155 sage: J1 = HadamardEJA(2)
2156 sage: J2 = RealSymmetricEJA(3)
2157 sage: J = DirectSumEJA(J1,J2)
2164 def __init__(self
, J1
, J2
, field
=AA
, **kwargs
):
2168 V
= VectorSpace(field
, n
)
2169 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2173 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2174 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2178 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2179 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2181 fdeja
= super(DirectSumEJA
, self
)
2182 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2183 self
.rank
.set_cache(J1
.rank() + J2
.rank())