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 # The multiplication table had better be square
102 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
107 self
.print_options(bracket
='')
109 # The multiplication table we're given is necessarily in terms
110 # of vectors, because we don't have an algebra yet for
111 # anything to be an element of. However, it's faster in the
112 # long run to have the multiplication table be in terms of
113 # algebra elements. We do this after calling the superclass
114 # constructor so that from_vector() knows what to do.
115 self
._multiplication
_table
= [
116 list(map(lambda x
: self
.from_vector(x
), ls
))
121 if not self
._is
_commutative
():
122 raise ValueError("algebra is not commutative")
123 if not self
._is
_jordanian
():
124 raise ValueError("Jordan identity does not hold")
125 if not self
._inner
_product
_is
_associative
():
126 raise ValueError("inner product is not associative")
128 def _element_constructor_(self
, elt
):
130 Construct an element of this algebra from its natural
133 This gets called only after the parent element _call_ method
134 fails to find a coercion for the argument.
138 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
140 ....: RealSymmetricEJA)
144 The identity in `S^n` is converted to the identity in the EJA::
146 sage: J = RealSymmetricEJA(3)
147 sage: I = matrix.identity(QQ,3)
148 sage: J(I) == J.one()
151 This skew-symmetric matrix can't be represented in the EJA::
153 sage: J = RealSymmetricEJA(3)
154 sage: A = matrix(QQ,3, lambda i,j: i-j)
156 Traceback (most recent call last):
158 ArithmeticError: vector is not in free module
162 Ensure that we can convert any element of the two non-matrix
163 simple algebras (whose natural representations are their usual
164 vector representations) back and forth faithfully::
166 sage: set_random_seed()
167 sage: J = HadamardEJA.random_instance()
168 sage: x = J.random_element()
169 sage: J(x.to_vector().column()) == x
171 sage: J = JordanSpinEJA.random_instance()
172 sage: x = J.random_element()
173 sage: J(x.to_vector().column()) == x
177 msg
= "not a naturally-represented algebra element"
179 # The superclass implementation of random_element()
180 # needs to be able to coerce "0" into the algebra.
182 elif elt
in self
.base_ring():
183 # Ensure that no base ring -> algebra coercion is performed
184 # by this method. There's some stupidity in sage that would
185 # otherwise propagate to this method; for example, sage thinks
186 # that the integer 3 belongs to the space of 2-by-2 matrices.
187 raise ValueError(msg
)
189 natural_basis
= self
.natural_basis()
190 basis_space
= natural_basis
[0].matrix_space()
191 if elt
not in basis_space
:
192 raise ValueError(msg
)
194 # Thanks for nothing! Matrix spaces aren't vector spaces in
195 # Sage, so we have to figure out its natural-basis coordinates
196 # ourselves. We use the basis space's ring instead of the
197 # element's ring because the basis space might be an algebraic
198 # closure whereas the base ring of the 3-by-3 identity matrix
199 # could be QQ instead of QQbar.
200 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
201 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
202 coords
= W
.coordinate_vector(_mat2vec(elt
))
203 return self
.from_vector(coords
)
206 def _max_test_case_size():
208 Return an integer "size" that is an upper bound on the size of
209 this algebra when it is used in a random test
210 case. Unfortunately, the term "size" is quite vague -- when
211 dealing with `R^n` under either the Hadamard or Jordan spin
212 product, the "size" refers to the dimension `n`. When dealing
213 with a matrix algebra (real symmetric or complex/quaternion
214 Hermitian), it refers to the size of the matrix, which is
215 far less than the dimension of the underlying vector space.
217 We default to five in this class, which is safe in `R^n`. The
218 matrix algebra subclasses (or any class where the "size" is
219 interpreted to be far less than the dimension) should override
220 with a smaller number.
226 Return a string representation of ``self``.
230 sage: from mjo.eja.eja_algebra import JordanSpinEJA
234 Ensure that it says what we think it says::
236 sage: JordanSpinEJA(2, field=AA)
237 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
238 sage: JordanSpinEJA(3, field=RDF)
239 Euclidean Jordan algebra of dimension 3 over Real Double Field
242 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
243 return fmt
.format(self
.dimension(), self
.base_ring())
245 def product_on_basis(self
, i
, j
):
246 return self
._multiplication
_table
[i
][j
]
248 def _is_commutative(self
):
250 Whether or not this algebra's multiplication table is commutative.
252 This method should of course always return ``True``, unless
253 this algebra was constructed with ``check=False`` and passed
254 an invalid multiplication table.
256 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
257 for i
in range(self
.dimension())
258 for j
in range(self
.dimension()) )
260 def _is_jordanian(self
):
262 Whether or not this algebra's multiplication table respects the
263 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
265 We only check one arrangement of `x` and `y`, so for a
266 ``True`` result to be truly true, you should also check
267 :meth:`_is_commutative`. This method should of course always
268 return ``True``, unless this algebra was constructed with
269 ``check=False`` and passed an invalid multiplication table.
271 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
273 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
274 for i
in range(self
.dimension())
275 for j
in range(self
.dimension()) )
277 def _inner_product_is_associative(self
):
279 Return whether or not this algebra's inner product `B` is
280 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
282 This method should of course always return ``True``, unless
283 this algebra was constructed with ``check=False`` and passed
284 an invalid multiplication table.
286 for i
in range(self
.dimension()):
287 for j
in range(self
.dimension()):
288 for k
in range(self
.dimension()):
292 if (x
*y
).inner_product(z
) != x
.inner_product(y
*z
):
298 def characteristic_polynomial_of(self
):
300 Return the algebra's "characteristic polynomial of" function,
301 which is itself a multivariate polynomial that, when evaluated
302 at the coordinates of some algebra element, returns that
303 element's characteristic polynomial.
305 The resulting polynomial has `n+1` variables, where `n` is the
306 dimension of this algebra. The first `n` variables correspond to
307 the coordinates of an algebra element: when evaluated at the
308 coordinates of an algebra element with respect to a certain
309 basis, the result is a univariate polynomial (in the one
310 remaining variable ``t``), namely the characteristic polynomial
315 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
319 The characteristic polynomial in the spin algebra is given in
320 Alizadeh, Example 11.11::
322 sage: J = JordanSpinEJA(3)
323 sage: p = J.characteristic_polynomial_of(); p
324 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
325 sage: xvec = J.one().to_vector()
329 By definition, the characteristic polynomial is a monic
330 degree-zero polynomial in a rank-zero algebra. Note that
331 Cayley-Hamilton is indeed satisfied since the polynomial
332 ``1`` evaluates to the identity element of the algebra on
335 sage: J = TrivialEJA()
336 sage: J.characteristic_polynomial_of()
343 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
344 a
= self
._charpoly
_coefficients
()
346 # We go to a bit of trouble here to reorder the
347 # indeterminates, so that it's easier to evaluate the
348 # characteristic polynomial at x's coordinates and get back
349 # something in terms of t, which is what we want.
350 S
= PolynomialRing(self
.base_ring(),'t')
354 S
= PolynomialRing(S
, R
.variable_names())
357 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
360 def inner_product(self
, x
, y
):
362 The inner product associated with this Euclidean Jordan algebra.
364 Defaults to the trace inner product, but can be overridden by
365 subclasses if they are sure that the necessary properties are
370 sage: from mjo.eja.eja_algebra import random_eja
374 Our inner product is "associative," which means the following for
375 a symmetric bilinear form::
377 sage: set_random_seed()
378 sage: J = random_eja()
379 sage: x,y,z = J.random_elements(3)
380 sage: (x*y).inner_product(z) == y.inner_product(x*z)
384 X
= x
.natural_representation()
385 Y
= y
.natural_representation()
386 return self
.natural_inner_product(X
,Y
)
389 def is_trivial(self
):
391 Return whether or not this algebra is trivial.
393 A trivial algebra contains only the zero element.
397 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
402 sage: J = ComplexHermitianEJA(3)
408 sage: J = TrivialEJA()
413 return self
.dimension() == 0
416 def multiplication_table(self
):
418 Return a visual representation of this algebra's multiplication
419 table (on basis elements).
423 sage: from mjo.eja.eja_algebra import JordanSpinEJA
427 sage: J = JordanSpinEJA(4)
428 sage: J.multiplication_table()
429 +----++----+----+----+----+
430 | * || e0 | e1 | e2 | e3 |
431 +====++====+====+====+====+
432 | e0 || e0 | e1 | e2 | e3 |
433 +----++----+----+----+----+
434 | e1 || e1 | e0 | 0 | 0 |
435 +----++----+----+----+----+
436 | e2 || e2 | 0 | e0 | 0 |
437 +----++----+----+----+----+
438 | e3 || e3 | 0 | 0 | e0 |
439 +----++----+----+----+----+
442 M
= list(self
._multiplication
_table
) # copy
443 for i
in range(len(M
)):
444 # M had better be "square"
445 M
[i
] = [self
.monomial(i
)] + M
[i
]
446 M
= [["*"] + list(self
.gens())] + M
447 return table(M
, header_row
=True, header_column
=True, frame
=True)
450 def natural_basis(self
):
452 Return a more-natural representation of this algebra's basis.
454 Every finite-dimensional Euclidean Jordan Algebra is a direct
455 sum of five simple algebras, four of which comprise Hermitian
456 matrices. This method returns the original "natural" basis
457 for our underlying vector space. (Typically, the natural basis
458 is used to construct the multiplication table in the first place.)
460 Note that this will always return a matrix. The standard basis
461 in `R^n` will be returned as `n`-by-`1` column matrices.
465 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
466 ....: RealSymmetricEJA)
470 sage: J = RealSymmetricEJA(2)
472 Finite family {0: e0, 1: e1, 2: e2}
473 sage: J.natural_basis()
475 [1 0] [ 0 0.7071067811865475?] [0 0]
476 [0 0], [0.7071067811865475? 0], [0 1]
481 sage: J = JordanSpinEJA(2)
483 Finite family {0: e0, 1: e1}
484 sage: J.natural_basis()
491 if self
._natural
_basis
is None:
492 M
= self
.natural_basis_space()
493 return tuple( M(b
.to_vector()) for b
in self
.basis() )
495 return self
._natural
_basis
498 def natural_basis_space(self
):
500 Return the matrix space in which this algebra's natural basis
503 Generally this will be an `n`-by-`1` column-vector space,
504 except when the algebra is trivial. There it's `n`-by-`n`
505 (where `n` is zero), to ensure that two elements of the
506 natural basis space (empty matrices) can be multiplied.
508 if self
.is_trivial():
509 return MatrixSpace(self
.base_ring(), 0)
510 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
511 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
513 return self
._natural
_basis
[0].matrix_space()
517 def natural_inner_product(X
,Y
):
519 Compute the inner product of two naturally-represented elements.
521 For example in the real symmetric matrix EJA, this will compute
522 the trace inner-product of two n-by-n symmetric matrices. The
523 default should work for the real cartesian product EJA, the
524 Jordan spin EJA, and the real symmetric matrices. The others
525 will have to be overridden.
527 return (X
.conjugate_transpose()*Y
).trace()
533 Return the unit element of this algebra.
537 sage: from mjo.eja.eja_algebra import (HadamardEJA,
542 sage: J = HadamardEJA(5)
544 e0 + e1 + e2 + e3 + e4
548 The identity element acts like the identity::
550 sage: set_random_seed()
551 sage: J = random_eja()
552 sage: x = J.random_element()
553 sage: J.one()*x == x and x*J.one() == x
556 The matrix of the unit element's operator is the identity::
558 sage: set_random_seed()
559 sage: J = random_eja()
560 sage: actual = J.one().operator().matrix()
561 sage: expected = matrix.identity(J.base_ring(), J.dimension())
562 sage: actual == expected
566 # We can brute-force compute the matrices of the operators
567 # that correspond to the basis elements of this algebra.
568 # If some linear combination of those basis elements is the
569 # algebra identity, then the same linear combination of
570 # their matrices has to be the identity matrix.
572 # Of course, matrices aren't vectors in sage, so we have to
573 # appeal to the "long vectors" isometry.
574 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
576 # Now we use basis linear algebra to find the coefficients,
577 # of the matrices-as-vectors-linear-combination, which should
578 # work for the original algebra basis too.
579 A
= matrix
.column(self
.base_ring(), oper_vecs
)
581 # We used the isometry on the left-hand side already, but we
582 # still need to do it for the right-hand side. Recall that we
583 # wanted something that summed to the identity matrix.
584 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
586 # Now if there's an identity element in the algebra, this should work.
587 coeffs
= A
.solve_right(b
)
588 return self
.linear_combination(zip(self
.gens(), coeffs
))
591 def peirce_decomposition(self
, c
):
593 The Peirce decomposition of this algebra relative to the
596 In the future, this can be extended to a complete system of
597 orthogonal idempotents.
601 - ``c`` -- an idempotent of this algebra.
605 A triple (J0, J5, J1) containing two subalgebras and one subspace
608 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
609 corresponding to the eigenvalue zero.
611 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
612 corresponding to the eigenvalue one-half.
614 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
615 corresponding to the eigenvalue one.
617 These are the only possible eigenspaces for that operator, and this
618 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
619 orthogonal, and are subalgebras of this algebra with the appropriate
624 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
628 The canonical example comes from the symmetric matrices, which
629 decompose into diagonal and off-diagonal parts::
631 sage: J = RealSymmetricEJA(3)
632 sage: C = matrix(QQ, [ [1,0,0],
636 sage: J0,J5,J1 = J.peirce_decomposition(c)
638 Euclidean Jordan algebra of dimension 1...
640 Vector space of degree 6 and dimension 2...
642 Euclidean Jordan algebra of dimension 3...
643 sage: J0.one().natural_representation()
647 sage: orig_df = AA.options.display_format
648 sage: AA.options.display_format = 'radical'
649 sage: J.from_vector(J5.basis()[0]).natural_representation()
653 sage: J.from_vector(J5.basis()[1]).natural_representation()
657 sage: AA.options.display_format = orig_df
658 sage: J1.one().natural_representation()
665 Every algebra decomposes trivially with respect to its identity
668 sage: set_random_seed()
669 sage: J = random_eja()
670 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
671 sage: J0.dimension() == 0 and J5.dimension() == 0
673 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
676 The decomposition is into eigenspaces, and its components are
677 therefore necessarily orthogonal. Moreover, the identity
678 elements in the two subalgebras are the projections onto their
679 respective subspaces of the superalgebra's identity element::
681 sage: set_random_seed()
682 sage: J = random_eja()
683 sage: x = J.random_element()
684 sage: if not J.is_trivial():
685 ....: while x.is_nilpotent():
686 ....: x = J.random_element()
687 sage: c = x.subalgebra_idempotent()
688 sage: J0,J5,J1 = J.peirce_decomposition(c)
690 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
691 ....: w = w.superalgebra_element()
692 ....: y = J.from_vector(y)
693 ....: z = z.superalgebra_element()
694 ....: ipsum += w.inner_product(y).abs()
695 ....: ipsum += w.inner_product(z).abs()
696 ....: ipsum += y.inner_product(z).abs()
699 sage: J1(c) == J1.one()
701 sage: J0(J.one() - c) == J0.one()
705 if not c
.is_idempotent():
706 raise ValueError("element is not idempotent: %s" % c
)
708 # Default these to what they should be if they turn out to be
709 # trivial, because eigenspaces_left() won't return eigenvalues
710 # corresponding to trivial spaces (e.g. it returns only the
711 # eigenspace corresponding to lambda=1 if you take the
712 # decomposition relative to the identity element).
713 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
714 J0
= trivial
# eigenvalue zero
715 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
716 J1
= trivial
# eigenvalue one
718 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
719 if eigval
== ~
(self
.base_ring()(2)):
722 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
723 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
729 raise ValueError("unexpected eigenvalue: %s" % eigval
)
734 def random_element(self
, thorough
=False):
736 Return a random element of this algebra.
738 Our algebra superclass method only returns a linear
739 combination of at most two basis elements. We instead
740 want the vector space "random element" method that
741 returns a more diverse selection.
745 - ``thorough`` -- (boolean; default False) whether or not we
746 should generate irrational coefficients for the random
747 element when our base ring is irrational; this slows the
748 algebra operations to a crawl, but any truly random method
752 # For a general base ring... maybe we can trust this to do the
753 # right thing? Unlikely, but.
754 V
= self
.vector_space()
755 v
= V
.random_element()
757 if self
.base_ring() is AA
:
758 # The "random element" method of the algebraic reals is
759 # stupid at the moment, and only returns integers between
760 # -2 and 2, inclusive:
762 # https://trac.sagemath.org/ticket/30875
764 # Instead, we implement our own "random vector" method,
765 # and then coerce that into the algebra. We use the vector
766 # space degree here instead of the dimension because a
767 # subalgebra could (for example) be spanned by only two
768 # vectors, each with five coordinates. We need to
769 # generate all five coordinates.
771 v
*= QQbar
.random_element().real()
773 v
*= QQ
.random_element()
775 return self
.from_vector(V
.coordinate_vector(v
))
777 def random_elements(self
, count
, thorough
=False):
779 Return ``count`` random elements as a tuple.
783 - ``thorough`` -- (boolean; default False) whether or not we
784 should generate irrational coefficients for the random
785 elements when our base ring is irrational; this slows the
786 algebra operations to a crawl, but any truly random method
791 sage: from mjo.eja.eja_algebra import JordanSpinEJA
795 sage: J = JordanSpinEJA(3)
796 sage: x,y,z = J.random_elements(3)
797 sage: all( [ x in J, y in J, z in J ])
799 sage: len( J.random_elements(10) ) == 10
803 return tuple( self
.random_element(thorough
)
804 for idx
in range(count
) )
807 def random_instance(cls
, field
=AA
, **kwargs
):
809 Return a random instance of this type of algebra.
811 Beware, this will crash for "most instances" because the
812 constructor below looks wrong.
814 if cls
is TrivialEJA
:
815 # The TrivialEJA class doesn't take an "n" argument because
819 n
= ZZ
.random_element(cls
._max
_test
_case
_size
() + 1)
820 return cls(n
, field
, **kwargs
)
823 def _charpoly_coefficients(self
):
825 The `r` polynomial coefficients of the "characteristic polynomial
829 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
830 R
= PolynomialRing(self
.base_ring(), var_names
)
832 F
= R
.fraction_field()
835 # From a result in my book, these are the entries of the
836 # basis representation of L_x.
837 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
840 L_x
= matrix(F
, n
, n
, L_x_i_j
)
843 if self
.rank
.is_in_cache():
845 # There's no need to pad the system with redundant
846 # columns if we *know* they'll be redundant.
849 # Compute an extra power in case the rank is equal to
850 # the dimension (otherwise, we would stop at x^(r-1)).
851 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
852 for k
in range(n
+1) ]
853 A
= matrix
.column(F
, x_powers
[:n
])
854 AE
= A
.extended_echelon_form()
861 # The theory says that only the first "r" coefficients are
862 # nonzero, and they actually live in the original polynomial
863 # ring and not the fraction field. We negate them because
864 # in the actual characteristic polynomial, they get moved
865 # to the other side where x^r lives.
866 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
871 Return the rank of this EJA.
873 This is a cached method because we know the rank a priori for
874 all of the algebras we can construct. Thus we can avoid the
875 expensive ``_charpoly_coefficients()`` call unless we truly
876 need to compute the whole characteristic polynomial.
880 sage: from mjo.eja.eja_algebra import (HadamardEJA,
882 ....: RealSymmetricEJA,
883 ....: ComplexHermitianEJA,
884 ....: QuaternionHermitianEJA,
889 The rank of the Jordan spin algebra is always two::
891 sage: JordanSpinEJA(2).rank()
893 sage: JordanSpinEJA(3).rank()
895 sage: JordanSpinEJA(4).rank()
898 The rank of the `n`-by-`n` Hermitian real, complex, or
899 quaternion matrices is `n`::
901 sage: RealSymmetricEJA(4).rank()
903 sage: ComplexHermitianEJA(3).rank()
905 sage: QuaternionHermitianEJA(2).rank()
910 Ensure that every EJA that we know how to construct has a
911 positive integer rank, unless the algebra is trivial in
912 which case its rank will be zero::
914 sage: set_random_seed()
915 sage: J = random_eja()
919 sage: r > 0 or (r == 0 and J.is_trivial())
922 Ensure that computing the rank actually works, since the ranks
923 of all simple algebras are known and will be cached by default::
925 sage: J = HadamardEJA(4)
926 sage: J.rank.clear_cache()
932 sage: J = JordanSpinEJA(4)
933 sage: J.rank.clear_cache()
939 sage: J = RealSymmetricEJA(3)
940 sage: J.rank.clear_cache()
946 sage: J = ComplexHermitianEJA(2)
947 sage: J.rank.clear_cache()
953 sage: J = QuaternionHermitianEJA(2)
954 sage: J.rank.clear_cache()
958 return len(self
._charpoly
_coefficients
())
961 def vector_space(self
):
963 Return the vector space that underlies this algebra.
967 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
971 sage: J = RealSymmetricEJA(2)
972 sage: J.vector_space()
973 Vector space of dimension 3 over...
976 return self
.zero().to_vector().parent().ambient_vector_space()
979 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
982 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
):
984 Return the Euclidean Jordan Algebra corresponding to the set
985 `R^n` under the Hadamard product.
987 Note: this is nothing more than the Cartesian product of ``n``
988 copies of the spin algebra. Once Cartesian product algebras
989 are implemented, this can go.
993 sage: from mjo.eja.eja_algebra import HadamardEJA
997 This multiplication table can be verified by hand::
999 sage: J = HadamardEJA(3)
1000 sage: e0,e1,e2 = J.gens()
1016 We can change the generator prefix::
1018 sage: HadamardEJA(3, prefix='r').gens()
1022 def __init__(self
, n
, field
=AA
, **kwargs
):
1023 V
= VectorSpace(field
, n
)
1024 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
1027 fdeja
= super(HadamardEJA
, self
)
1028 fdeja
.__init
__(field
, mult_table
, **kwargs
)
1029 self
.rank
.set_cache(n
)
1031 def inner_product(self
, x
, y
):
1033 Faster to reimplement than to use natural representations.
1037 sage: from mjo.eja.eja_algebra import HadamardEJA
1041 Ensure that this is the usual inner product for the algebras
1044 sage: set_random_seed()
1045 sage: J = HadamardEJA.random_instance()
1046 sage: x,y = J.random_elements(2)
1047 sage: X = x.natural_representation()
1048 sage: Y = y.natural_representation()
1049 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1053 return x
.to_vector().inner_product(y
.to_vector())
1056 def random_eja(field
=AA
):
1058 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1062 sage: from mjo.eja.eja_algebra import random_eja
1067 Euclidean Jordan algebra of dimension...
1070 classname
= choice([TrivialEJA
,
1074 ComplexHermitianEJA
,
1075 QuaternionHermitianEJA
])
1076 return classname
.random_instance(field
=field
)
1081 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1083 def _max_test_case_size():
1084 # Play it safe, since this will be squared and the underlying
1085 # field can have dimension 4 (quaternions) too.
1088 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1090 Compared to the superclass constructor, we take a basis instead of
1091 a multiplication table because the latter can be computed in terms
1092 of the former when the product is known (like it is here).
1094 # Used in this class's fast _charpoly_coefficients() override.
1095 self
._basis
_normalizers
= None
1097 # We're going to loop through this a few times, so now's a good
1098 # time to ensure that it isn't a generator expression.
1099 basis
= tuple(basis
)
1101 if len(basis
) > 1 and normalize_basis
:
1102 # We'll need sqrt(2) to normalize the basis, and this
1103 # winds up in the multiplication table, so the whole
1104 # algebra needs to be over the field extension.
1105 R
= PolynomialRing(field
, 'z')
1108 if p
.is_irreducible():
1109 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1110 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1111 self
._basis
_normalizers
= tuple(
1112 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1113 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1115 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1117 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
1118 fdeja
.__init
__(field
, Qs
, natural_basis
=basis
, **kwargs
)
1123 def _charpoly_coefficients(self
):
1125 Override the parent method with something that tries to compute
1126 over a faster (non-extension) field.
1128 if self
._basis
_normalizers
is None:
1129 # We didn't normalize, so assume that the basis we started
1130 # with had entries in a nice field.
1131 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1133 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1134 self
._basis
_normalizers
) )
1136 # Do this over the rationals and convert back at the end.
1137 # Only works because we know the entries of the basis are
1139 J
= MatrixEuclideanJordanAlgebra(QQ
,
1141 normalize_basis
=False)
1142 a
= J
._charpoly
_coefficients
()
1144 # Unfortunately, changing the basis does change the
1145 # coefficients of the characteristic polynomial, but since
1146 # these are really the coefficients of the "characteristic
1147 # polynomial of" function, everything is still nice and
1148 # unevaluated. It's therefore "obvious" how scaling the
1149 # basis affects the coordinate variables X1, X2, et
1150 # cetera. Scaling the first basis vector up by "n" adds a
1151 # factor of 1/n into every "X1" term, for example. So here
1152 # we simply undo the basis_normalizer scaling that we
1153 # performed earlier.
1155 # The a[0] access here is safe because trivial algebras
1156 # won't have any basis normalizers and therefore won't
1157 # make it to this "else" branch.
1158 XS
= a
[0].parent().gens()
1159 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1160 for i
in range(len(XS
)) }
1161 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1165 def multiplication_table_from_matrix_basis(basis
):
1167 At least three of the five simple Euclidean Jordan algebras have the
1168 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1169 multiplication on the right is matrix multiplication. Given a basis
1170 for the underlying matrix space, this function returns a
1171 multiplication table (obtained by looping through the basis
1172 elements) for an algebra of those matrices.
1174 # In S^2, for example, we nominally have four coordinates even
1175 # though the space is of dimension three only. The vector space V
1176 # is supposed to hold the entire long vector, and the subspace W
1177 # of V will be spanned by the vectors that arise from symmetric
1178 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1182 field
= basis
[0].base_ring()
1183 dimension
= basis
[0].nrows()
1185 V
= VectorSpace(field
, dimension
**2)
1186 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1188 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1191 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1192 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1200 Embed the matrix ``M`` into a space of real matrices.
1202 The matrix ``M`` can have entries in any field at the moment:
1203 the real numbers, complex numbers, or quaternions. And although
1204 they are not a field, we can probably support octonions at some
1205 point, too. This function returns a real matrix that "acts like"
1206 the original with respect to matrix multiplication; i.e.
1208 real_embed(M*N) = real_embed(M)*real_embed(N)
1211 raise NotImplementedError
1215 def real_unembed(M
):
1217 The inverse of :meth:`real_embed`.
1219 raise NotImplementedError
1223 def natural_inner_product(cls
,X
,Y
):
1224 Xu
= cls
.real_unembed(X
)
1225 Yu
= cls
.real_unembed(Y
)
1226 tr
= (Xu
*Yu
).trace()
1229 # It's real already.
1232 # Otherwise, try the thing that works for complex numbers; and
1233 # if that doesn't work, the thing that works for quaternions.
1235 return tr
.vector()[0] # real part, imag part is index 1
1236 except AttributeError:
1237 # A quaternions doesn't have a vector() method, but does
1238 # have coefficient_tuple() method that returns the
1239 # coefficients of 1, i, j, and k -- in that order.
1240 return tr
.coefficient_tuple()[0]
1243 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1247 The identity function, for embedding real matrices into real
1253 def real_unembed(M
):
1255 The identity function, for unembedding real matrices from real
1261 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1263 The rank-n simple EJA consisting of real symmetric n-by-n
1264 matrices, the usual symmetric Jordan product, and the trace inner
1265 product. It has dimension `(n^2 + n)/2` over the reals.
1269 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1273 sage: J = RealSymmetricEJA(2)
1274 sage: e0, e1, e2 = J.gens()
1282 In theory, our "field" can be any subfield of the reals::
1284 sage: RealSymmetricEJA(2, RDF)
1285 Euclidean Jordan algebra of dimension 3 over Real Double Field
1286 sage: RealSymmetricEJA(2, RR)
1287 Euclidean Jordan algebra of dimension 3 over Real Field with
1288 53 bits of precision
1292 The dimension of this algebra is `(n^2 + n) / 2`::
1294 sage: set_random_seed()
1295 sage: n_max = RealSymmetricEJA._max_test_case_size()
1296 sage: n = ZZ.random_element(1, n_max)
1297 sage: J = RealSymmetricEJA(n)
1298 sage: J.dimension() == (n^2 + n)/2
1301 The Jordan multiplication is what we think it is::
1303 sage: set_random_seed()
1304 sage: J = RealSymmetricEJA.random_instance()
1305 sage: x,y = J.random_elements(2)
1306 sage: actual = (x*y).natural_representation()
1307 sage: X = x.natural_representation()
1308 sage: Y = y.natural_representation()
1309 sage: expected = (X*Y + Y*X)/2
1310 sage: actual == expected
1312 sage: J(expected) == x*y
1315 We can change the generator prefix::
1317 sage: RealSymmetricEJA(3, prefix='q').gens()
1318 (q0, q1, q2, q3, q4, q5)
1320 Our natural basis is normalized with respect to the natural inner
1321 product unless we specify otherwise::
1323 sage: set_random_seed()
1324 sage: J = RealSymmetricEJA.random_instance()
1325 sage: all( b.norm() == 1 for b in J.gens() )
1328 Since our natural basis is normalized with respect to the natural
1329 inner product, and since we know that this algebra is an EJA, any
1330 left-multiplication operator's matrix will be symmetric because
1331 natural->EJA basis representation is an isometry and within the EJA
1332 the operator is self-adjoint by the Jordan axiom::
1334 sage: set_random_seed()
1335 sage: x = RealSymmetricEJA.random_instance().random_element()
1336 sage: x.operator().matrix().is_symmetric()
1339 We can construct the (trivial) algebra of rank zero::
1341 sage: RealSymmetricEJA(0)
1342 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1346 def _denormalized_basis(cls
, n
, field
):
1348 Return a basis for the space of real symmetric n-by-n matrices.
1352 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1356 sage: set_random_seed()
1357 sage: n = ZZ.random_element(1,5)
1358 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1359 sage: all( M.is_symmetric() for M in B)
1363 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1367 for j
in range(i
+1):
1368 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1372 Sij
= Eij
+ Eij
.transpose()
1378 def _max_test_case_size():
1379 return 4 # Dimension 10
1382 def __init__(self
, n
, field
=AA
, **kwargs
):
1383 basis
= self
._denormalized
_basis
(n
, field
)
1384 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, **kwargs
)
1385 self
.rank
.set_cache(n
)
1388 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1392 Embed the n-by-n complex matrix ``M`` into the space of real
1393 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1394 bi` to the block matrix ``[[a,b],[-b,a]]``.
1398 sage: from mjo.eja.eja_algebra import \
1399 ....: ComplexMatrixEuclideanJordanAlgebra
1403 sage: F = QuadraticField(-1, 'I')
1404 sage: x1 = F(4 - 2*i)
1405 sage: x2 = F(1 + 2*i)
1408 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1409 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1418 Embedding is a homomorphism (isomorphism, in fact)::
1420 sage: set_random_seed()
1421 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1422 sage: n = ZZ.random_element(n_max)
1423 sage: F = QuadraticField(-1, 'I')
1424 sage: X = random_matrix(F, n)
1425 sage: Y = random_matrix(F, n)
1426 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1427 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1428 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1435 raise ValueError("the matrix 'M' must be square")
1437 # We don't need any adjoined elements...
1438 field
= M
.base_ring().base_ring()
1442 a
= z
.list()[0] # real part, I guess
1443 b
= z
.list()[1] # imag part, I guess
1444 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1446 return matrix
.block(field
, n
, blocks
)
1450 def real_unembed(M
):
1452 The inverse of _embed_complex_matrix().
1456 sage: from mjo.eja.eja_algebra import \
1457 ....: ComplexMatrixEuclideanJordanAlgebra
1461 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1462 ....: [-2, 1, -4, 3],
1463 ....: [ 9, 10, 11, 12],
1464 ....: [-10, 9, -12, 11] ])
1465 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1467 [ 10*I + 9 12*I + 11]
1471 Unembedding is the inverse of embedding::
1473 sage: set_random_seed()
1474 sage: F = QuadraticField(-1, 'I')
1475 sage: M = random_matrix(F, 3)
1476 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1477 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1483 raise ValueError("the matrix 'M' must be square")
1484 if not n
.mod(2).is_zero():
1485 raise ValueError("the matrix 'M' must be a complex embedding")
1487 # If "M" was normalized, its base ring might have roots
1488 # adjoined and they can stick around after unembedding.
1489 field
= M
.base_ring()
1490 R
= PolynomialRing(field
, 'z')
1493 # Sage doesn't know how to embed AA into QQbar, i.e. how
1494 # to adjoin sqrt(-1) to AA.
1497 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1500 # Go top-left to bottom-right (reading order), converting every
1501 # 2-by-2 block we see to a single complex element.
1503 for k
in range(n
/2):
1504 for j
in range(n
/2):
1505 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1506 if submat
[0,0] != submat
[1,1]:
1507 raise ValueError('bad on-diagonal submatrix')
1508 if submat
[0,1] != -submat
[1,0]:
1509 raise ValueError('bad off-diagonal submatrix')
1510 z
= submat
[0,0] + submat
[0,1]*i
1513 return matrix(F
, n
/2, elements
)
1517 def natural_inner_product(cls
,X
,Y
):
1519 Compute a natural inner product in this algebra directly from
1524 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1528 This gives the same answer as the slow, default method implemented
1529 in :class:`MatrixEuclideanJordanAlgebra`::
1531 sage: set_random_seed()
1532 sage: J = ComplexHermitianEJA.random_instance()
1533 sage: x,y = J.random_elements(2)
1534 sage: Xe = x.natural_representation()
1535 sage: Ye = y.natural_representation()
1536 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1537 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1538 sage: expected = (X*Y).trace().real()
1539 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1540 sage: actual == expected
1544 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1547 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1549 The rank-n simple EJA consisting of complex Hermitian n-by-n
1550 matrices over the real numbers, the usual symmetric Jordan product,
1551 and the real-part-of-trace inner product. It has dimension `n^2` over
1556 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1560 In theory, our "field" can be any subfield of the reals::
1562 sage: ComplexHermitianEJA(2, RDF)
1563 Euclidean Jordan algebra of dimension 4 over Real Double Field
1564 sage: ComplexHermitianEJA(2, RR)
1565 Euclidean Jordan algebra of dimension 4 over Real Field with
1566 53 bits of precision
1570 The dimension of this algebra is `n^2`::
1572 sage: set_random_seed()
1573 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1574 sage: n = ZZ.random_element(1, n_max)
1575 sage: J = ComplexHermitianEJA(n)
1576 sage: J.dimension() == n^2
1579 The Jordan multiplication is what we think it is::
1581 sage: set_random_seed()
1582 sage: J = ComplexHermitianEJA.random_instance()
1583 sage: x,y = J.random_elements(2)
1584 sage: actual = (x*y).natural_representation()
1585 sage: X = x.natural_representation()
1586 sage: Y = y.natural_representation()
1587 sage: expected = (X*Y + Y*X)/2
1588 sage: actual == expected
1590 sage: J(expected) == x*y
1593 We can change the generator prefix::
1595 sage: ComplexHermitianEJA(2, prefix='z').gens()
1598 Our natural basis is normalized with respect to the natural inner
1599 product unless we specify otherwise::
1601 sage: set_random_seed()
1602 sage: J = ComplexHermitianEJA.random_instance()
1603 sage: all( b.norm() == 1 for b in J.gens() )
1606 Since our natural basis is normalized with respect to the natural
1607 inner product, and since we know that this algebra is an EJA, any
1608 left-multiplication operator's matrix will be symmetric because
1609 natural->EJA basis representation is an isometry and within the EJA
1610 the operator is self-adjoint by the Jordan axiom::
1612 sage: set_random_seed()
1613 sage: x = ComplexHermitianEJA.random_instance().random_element()
1614 sage: x.operator().matrix().is_symmetric()
1617 We can construct the (trivial) algebra of rank zero::
1619 sage: ComplexHermitianEJA(0)
1620 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1625 def _denormalized_basis(cls
, n
, field
):
1627 Returns a basis for the space of complex Hermitian n-by-n matrices.
1629 Why do we embed these? Basically, because all of numerical linear
1630 algebra assumes that you're working with vectors consisting of `n`
1631 entries from a field and scalars from the same field. There's no way
1632 to tell SageMath that (for example) the vectors contain complex
1633 numbers, while the scalar field is real.
1637 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1641 sage: set_random_seed()
1642 sage: n = ZZ.random_element(1,5)
1643 sage: field = QuadraticField(2, 'sqrt2')
1644 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1645 sage: all( M.is_symmetric() for M in B)
1649 R
= PolynomialRing(field
, 'z')
1651 F
= field
.extension(z
**2 + 1, 'I')
1654 # This is like the symmetric case, but we need to be careful:
1656 # * We want conjugate-symmetry, not just symmetry.
1657 # * The diagonal will (as a result) be real.
1661 for j
in range(i
+1):
1662 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1664 Sij
= cls
.real_embed(Eij
)
1667 # The second one has a minus because it's conjugated.
1668 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1670 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1673 # Since we embedded these, we can drop back to the "field" that we
1674 # started with instead of the complex extension "F".
1675 return ( s
.change_ring(field
) for s
in S
)
1678 def __init__(self
, n
, field
=AA
, **kwargs
):
1679 basis
= self
._denormalized
_basis
(n
,field
)
1680 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1681 self
.rank
.set_cache(n
)
1684 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1688 Embed the n-by-n quaternion matrix ``M`` into the space of real
1689 matrices of size 4n-by-4n by first sending each quaternion entry `z
1690 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1691 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1696 sage: from mjo.eja.eja_algebra import \
1697 ....: QuaternionMatrixEuclideanJordanAlgebra
1701 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1702 sage: i,j,k = Q.gens()
1703 sage: x = 1 + 2*i + 3*j + 4*k
1704 sage: M = matrix(Q, 1, [[x]])
1705 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1711 Embedding is a homomorphism (isomorphism, in fact)::
1713 sage: set_random_seed()
1714 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1715 sage: n = ZZ.random_element(n_max)
1716 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1717 sage: X = random_matrix(Q, n)
1718 sage: Y = random_matrix(Q, n)
1719 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1720 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1721 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1726 quaternions
= M
.base_ring()
1729 raise ValueError("the matrix 'M' must be square")
1731 F
= QuadraticField(-1, 'I')
1736 t
= z
.coefficient_tuple()
1741 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1742 [-c
+ d
*i
, a
- b
*i
]])
1743 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1744 blocks
.append(realM
)
1746 # We should have real entries by now, so use the realest field
1747 # we've got for the return value.
1748 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1753 def real_unembed(M
):
1755 The inverse of _embed_quaternion_matrix().
1759 sage: from mjo.eja.eja_algebra import \
1760 ....: QuaternionMatrixEuclideanJordanAlgebra
1764 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1765 ....: [-2, 1, -4, 3],
1766 ....: [-3, 4, 1, -2],
1767 ....: [-4, -3, 2, 1]])
1768 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1769 [1 + 2*i + 3*j + 4*k]
1773 Unembedding is the inverse of embedding::
1775 sage: set_random_seed()
1776 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1777 sage: M = random_matrix(Q, 3)
1778 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1779 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1785 raise ValueError("the matrix 'M' must be square")
1786 if not n
.mod(4).is_zero():
1787 raise ValueError("the matrix 'M' must be a quaternion embedding")
1789 # Use the base ring of the matrix to ensure that its entries can be
1790 # multiplied by elements of the quaternion algebra.
1791 field
= M
.base_ring()
1792 Q
= QuaternionAlgebra(field
,-1,-1)
1795 # Go top-left to bottom-right (reading order), converting every
1796 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1799 for l
in range(n
/4):
1800 for m
in range(n
/4):
1801 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1802 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1803 if submat
[0,0] != submat
[1,1].conjugate():
1804 raise ValueError('bad on-diagonal submatrix')
1805 if submat
[0,1] != -submat
[1,0].conjugate():
1806 raise ValueError('bad off-diagonal submatrix')
1807 z
= submat
[0,0].real()
1808 z
+= submat
[0,0].imag()*i
1809 z
+= submat
[0,1].real()*j
1810 z
+= submat
[0,1].imag()*k
1813 return matrix(Q
, n
/4, elements
)
1817 def natural_inner_product(cls
,X
,Y
):
1819 Compute a natural inner product in this algebra directly from
1824 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1828 This gives the same answer as the slow, default method implemented
1829 in :class:`MatrixEuclideanJordanAlgebra`::
1831 sage: set_random_seed()
1832 sage: J = QuaternionHermitianEJA.random_instance()
1833 sage: x,y = J.random_elements(2)
1834 sage: Xe = x.natural_representation()
1835 sage: Ye = y.natural_representation()
1836 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1837 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1838 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1839 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1840 sage: actual == expected
1844 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1847 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1849 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1850 matrices, the usual symmetric Jordan product, and the
1851 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1856 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1860 In theory, our "field" can be any subfield of the reals::
1862 sage: QuaternionHermitianEJA(2, RDF)
1863 Euclidean Jordan algebra of dimension 6 over Real Double Field
1864 sage: QuaternionHermitianEJA(2, RR)
1865 Euclidean Jordan algebra of dimension 6 over Real Field with
1866 53 bits of precision
1870 The dimension of this algebra is `2*n^2 - n`::
1872 sage: set_random_seed()
1873 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1874 sage: n = ZZ.random_element(1, n_max)
1875 sage: J = QuaternionHermitianEJA(n)
1876 sage: J.dimension() == 2*(n^2) - n
1879 The Jordan multiplication is what we think it is::
1881 sage: set_random_seed()
1882 sage: J = QuaternionHermitianEJA.random_instance()
1883 sage: x,y = J.random_elements(2)
1884 sage: actual = (x*y).natural_representation()
1885 sage: X = x.natural_representation()
1886 sage: Y = y.natural_representation()
1887 sage: expected = (X*Y + Y*X)/2
1888 sage: actual == expected
1890 sage: J(expected) == x*y
1893 We can change the generator prefix::
1895 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1896 (a0, a1, a2, a3, a4, a5)
1898 Our natural basis is normalized with respect to the natural inner
1899 product unless we specify otherwise::
1901 sage: set_random_seed()
1902 sage: J = QuaternionHermitianEJA.random_instance()
1903 sage: all( b.norm() == 1 for b in J.gens() )
1906 Since our natural basis is normalized with respect to the natural
1907 inner product, and since we know that this algebra is an EJA, any
1908 left-multiplication operator's matrix will be symmetric because
1909 natural->EJA basis representation is an isometry and within the EJA
1910 the operator is self-adjoint by the Jordan axiom::
1912 sage: set_random_seed()
1913 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1914 sage: x.operator().matrix().is_symmetric()
1917 We can construct the (trivial) algebra of rank zero::
1919 sage: QuaternionHermitianEJA(0)
1920 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1924 def _denormalized_basis(cls
, n
, field
):
1926 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1928 Why do we embed these? Basically, because all of numerical
1929 linear algebra assumes that you're working with vectors consisting
1930 of `n` entries from a field and scalars from the same field. There's
1931 no way to tell SageMath that (for example) the vectors contain
1932 complex numbers, while the scalar field is real.
1936 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1940 sage: set_random_seed()
1941 sage: n = ZZ.random_element(1,5)
1942 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1943 sage: all( M.is_symmetric() for M in B )
1947 Q
= QuaternionAlgebra(QQ
,-1,-1)
1950 # This is like the symmetric case, but we need to be careful:
1952 # * We want conjugate-symmetry, not just symmetry.
1953 # * The diagonal will (as a result) be real.
1957 for j
in range(i
+1):
1958 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1960 Sij
= cls
.real_embed(Eij
)
1963 # The second, third, and fourth ones have a minus
1964 # because they're conjugated.
1965 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1967 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1969 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1971 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1974 # Since we embedded these, we can drop back to the "field" that we
1975 # started with instead of the quaternion algebra "Q".
1976 return ( s
.change_ring(field
) for s
in S
)
1979 def __init__(self
, n
, field
=AA
, **kwargs
):
1980 basis
= self
._denormalized
_basis
(n
,field
)
1981 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1982 self
.rank
.set_cache(n
)
1985 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1987 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1988 with the half-trace inner product and jordan product ``x*y =
1989 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
1990 symmetric positive-definite "bilinear form" matrix. It has
1991 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
1992 when ``B`` is the identity matrix of order ``n-1``.
1996 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1997 ....: JordanSpinEJA)
2001 When no bilinear form is specified, the identity matrix is used,
2002 and the resulting algebra is the Jordan spin algebra::
2004 sage: J0 = BilinearFormEJA(3)
2005 sage: J1 = JordanSpinEJA(3)
2006 sage: J0.multiplication_table() == J0.multiplication_table()
2011 We can create a zero-dimensional algebra::
2013 sage: J = BilinearFormEJA(0)
2017 We can check the multiplication condition given in the Jordan, von
2018 Neumann, and Wigner paper (and also discussed on my "On the
2019 symmetry..." paper). Note that this relies heavily on the standard
2020 choice of basis, as does anything utilizing the bilinear form matrix::
2022 sage: set_random_seed()
2023 sage: n = ZZ.random_element(5)
2024 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2025 sage: B = M.transpose()*M
2026 sage: J = BilinearFormEJA(n, B=B)
2027 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2028 sage: V = J.vector_space()
2029 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2030 ....: for ei in eis ]
2031 sage: actual = [ sis[i]*sis[j]
2032 ....: for i in range(n-1)
2033 ....: for j in range(n-1) ]
2034 sage: expected = [ J.one() if i == j else J.zero()
2035 ....: for i in range(n-1)
2036 ....: for j in range(n-1) ]
2037 sage: actual == expected
2040 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
2042 self
._B
= matrix
.identity(field
, max(0,n
-1))
2046 V
= VectorSpace(field
, n
)
2047 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2056 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
2057 zbar
= y0
*xbar
+ x0
*ybar
2058 z
= V([z0
] + zbar
.list())
2059 mult_table
[i
][j
] = z
2061 # The rank of this algebra is two, unless we're in a
2062 # one-dimensional ambient space (because the rank is bounded
2063 # by the ambient dimension).
2064 fdeja
= super(BilinearFormEJA
, self
)
2065 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2066 self
.rank
.set_cache(min(n
,2))
2068 def inner_product(self
, x
, y
):
2070 Half of the trace inner product.
2072 This is defined so that the special case of the Jordan spin
2073 algebra gets the usual inner product.
2077 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2081 Ensure that this is one-half of the trace inner-product when
2082 the algebra isn't just the reals (when ``n`` isn't one). This
2083 is in Faraut and Koranyi, and also my "On the symmetry..."
2086 sage: set_random_seed()
2087 sage: n = ZZ.random_element(2,5)
2088 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2089 sage: B = M.transpose()*M
2090 sage: J = BilinearFormEJA(n, B=B)
2091 sage: x = J.random_element()
2092 sage: y = J.random_element()
2093 sage: x.inner_product(y) == (x*y).trace()/2
2097 xvec
= x
.to_vector()
2099 yvec
= y
.to_vector()
2101 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
2104 class JordanSpinEJA(BilinearFormEJA
):
2106 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2107 with the usual inner product and jordan product ``x*y =
2108 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2113 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2117 This multiplication table can be verified by hand::
2119 sage: J = JordanSpinEJA(4)
2120 sage: e0,e1,e2,e3 = J.gens()
2136 We can change the generator prefix::
2138 sage: JordanSpinEJA(2, prefix='B').gens()
2143 Ensure that we have the usual inner product on `R^n`::
2145 sage: set_random_seed()
2146 sage: J = JordanSpinEJA.random_instance()
2147 sage: x,y = J.random_elements(2)
2148 sage: X = x.natural_representation()
2149 sage: Y = y.natural_representation()
2150 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2154 def __init__(self
, n
, field
=AA
, **kwargs
):
2155 # This is a special case of the BilinearFormEJA with the identity
2156 # matrix as its bilinear form.
2157 return super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2160 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2162 The trivial Euclidean Jordan algebra consisting of only a zero element.
2166 sage: from mjo.eja.eja_algebra import TrivialEJA
2170 sage: J = TrivialEJA()
2177 sage: 7*J.one()*12*J.one()
2179 sage: J.one().inner_product(J.one())
2181 sage: J.one().norm()
2183 sage: J.one().subalgebra_generated_by()
2184 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2189 def __init__(self
, field
=AA
, **kwargs
):
2191 fdeja
= super(TrivialEJA
, self
)
2192 # The rank is zero using my definition, namely the dimension of the
2193 # largest subalgebra generated by any element.
2194 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2195 self
.rank
.set_cache(0)
2198 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2200 The external (orthogonal) direct sum of two other Euclidean Jordan
2201 algebras. Essentially the Cartesian product of its two factors.
2202 Every Euclidean Jordan algebra decomposes into an orthogonal
2203 direct sum of simple Euclidean Jordan algebras, so no generality
2204 is lost by providing only this construction.
2208 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2209 ....: RealSymmetricEJA,
2214 sage: J1 = HadamardEJA(2)
2215 sage: J2 = RealSymmetricEJA(3)
2216 sage: J = DirectSumEJA(J1,J2)
2223 def __init__(self
, J1
, J2
, field
=AA
, **kwargs
):
2227 V
= VectorSpace(field
, n
)
2228 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2232 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2233 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2237 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2238 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2240 fdeja
= super(DirectSumEJA
, self
)
2241 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2242 self
.rank
.set_cache(J1
.rank() + J2
.rank())