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
]
238 def _a_regular_element(self
):
240 Guess a regular element. Needed to compute the basis for our
241 characteristic polynomial coefficients.
245 sage: from mjo.eja.eja_algebra import random_eja
249 Ensure that this hacky method succeeds for every algebra that we
250 know how to construct::
252 sage: set_random_seed()
253 sage: J = random_eja()
254 sage: J._a_regular_element().is_regular()
259 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
260 if not z
.is_regular():
261 raise ValueError("don't know a regular element")
266 def _charpoly_basis_space(self
):
268 Return the vector space spanned by the basis used in our
269 characteristic polynomial coefficients. This is used not only to
270 compute those coefficients, but also any time we need to
271 evaluate the coefficients (like when we compute the trace or
274 z
= self
._a
_regular
_element
()
275 # Don't use the parent vector space directly here in case this
276 # happens to be a subalgebra. In that case, we would be e.g.
277 # two-dimensional but span_of_basis() would expect three
279 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
280 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
281 V1
= V
.span_of_basis( basis
)
282 b
= (V1
.basis() + V1
.complement().basis())
283 return V
.span_of_basis(b
)
286 def _charpoly_coeff(self
, i
):
288 Return the coefficient polynomial "a_{i}" of this algebra's
289 general characteristic polynomial.
291 Having this be a separate cached method lets us compute and
292 store the trace/determinant (a_{r-1} and a_{0} respectively)
293 separate from the entire characteristic polynomial.
295 return self
._charpoly
_coefficients
()[i
]
299 def characteristic_polynomial(self
):
301 Return a characteristic polynomial that works for all elements
304 The resulting polynomial has `n+1` variables, where `n` is the
305 dimension of this algebra. The first `n` variables correspond to
306 the coordinates of an algebra element: when evaluated at the
307 coordinates of an algebra element with respect to a certain
308 basis, the result is a univariate polynomial (in the one
309 remaining variable ``t``), namely the characteristic polynomial
314 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
318 The characteristic polynomial in the spin algebra is given in
319 Alizadeh, Example 11.11::
321 sage: J = JordanSpinEJA(3)
322 sage: p = J.characteristic_polynomial(); p
323 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
324 sage: xvec = J.one().to_vector()
328 By definition, the characteristic polynomial is a monic
329 degree-zero polynomial in a rank-zero algebra. Note that
330 Cayley-Hamilton is indeed satisfied since the polynomial
331 ``1`` evaluates to the identity element of the algebra on
334 sage: J = TrivialEJA()
335 sage: J.characteristic_polynomial()
342 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
343 a
= self
._charpoly
_coefficients
()
345 # We go to a bit of trouble here to reorder the
346 # indeterminates, so that it's easier to evaluate the
347 # characteristic polynomial at x's coordinates and get back
348 # something in terms of t, which is what we want.
349 S
= PolynomialRing(self
.base_ring(),'t')
353 S
= PolynomialRing(S
, R
.variable_names())
356 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
359 def inner_product(self
, x
, y
):
361 The inner product associated with this Euclidean Jordan algebra.
363 Defaults to the trace inner product, but can be overridden by
364 subclasses if they are sure that the necessary properties are
369 sage: from mjo.eja.eja_algebra import random_eja
373 Our inner product is "associative," which means the following for
374 a symmetric bilinear form::
376 sage: set_random_seed()
377 sage: J = random_eja()
378 sage: x,y,z = J.random_elements(3)
379 sage: (x*y).inner_product(z) == y.inner_product(x*z)
383 X
= x
.natural_representation()
384 Y
= y
.natural_representation()
385 return self
.natural_inner_product(X
,Y
)
388 def is_trivial(self
):
390 Return whether or not this algebra is trivial.
392 A trivial algebra contains only the zero element.
396 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
401 sage: J = ComplexHermitianEJA(3)
407 sage: J = TrivialEJA()
412 return self
.dimension() == 0
415 def multiplication_table(self
):
417 Return a visual representation of this algebra's multiplication
418 table (on basis elements).
422 sage: from mjo.eja.eja_algebra import JordanSpinEJA
426 sage: J = JordanSpinEJA(4)
427 sage: J.multiplication_table()
428 +----++----+----+----+----+
429 | * || e0 | e1 | e2 | e3 |
430 +====++====+====+====+====+
431 | e0 || e0 | e1 | e2 | e3 |
432 +----++----+----+----+----+
433 | e1 || e1 | e0 | 0 | 0 |
434 +----++----+----+----+----+
435 | e2 || e2 | 0 | e0 | 0 |
436 +----++----+----+----+----+
437 | e3 || e3 | 0 | 0 | e0 |
438 +----++----+----+----+----+
441 M
= list(self
._multiplication
_table
) # copy
442 for i
in range(len(M
)):
443 # M had better be "square"
444 M
[i
] = [self
.monomial(i
)] + M
[i
]
445 M
= [["*"] + list(self
.gens())] + M
446 return table(M
, header_row
=True, header_column
=True, frame
=True)
449 def natural_basis(self
):
451 Return a more-natural representation of this algebra's basis.
453 Every finite-dimensional Euclidean Jordan Algebra is a direct
454 sum of five simple algebras, four of which comprise Hermitian
455 matrices. This method returns the original "natural" basis
456 for our underlying vector space. (Typically, the natural basis
457 is used to construct the multiplication table in the first place.)
459 Note that this will always return a matrix. The standard basis
460 in `R^n` will be returned as `n`-by-`1` column matrices.
464 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
465 ....: RealSymmetricEJA)
469 sage: J = RealSymmetricEJA(2)
471 Finite family {0: e0, 1: e1, 2: e2}
472 sage: J.natural_basis()
474 [1 0] [ 0 0.7071067811865475?] [0 0]
475 [0 0], [0.7071067811865475? 0], [0 1]
480 sage: J = JordanSpinEJA(2)
482 Finite family {0: e0, 1: e1}
483 sage: J.natural_basis()
490 if self
._natural
_basis
is None:
491 M
= self
.natural_basis_space()
492 return tuple( M(b
.to_vector()) for b
in self
.basis() )
494 return self
._natural
_basis
497 def natural_basis_space(self
):
499 Return the matrix space in which this algebra's natural basis
502 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
503 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
505 return self
._natural
_basis
[0].matrix_space()
509 def natural_inner_product(X
,Y
):
511 Compute the inner product of two naturally-represented elements.
513 For example in the real symmetric matrix EJA, this will compute
514 the trace inner-product of two n-by-n symmetric matrices. The
515 default should work for the real cartesian product EJA, the
516 Jordan spin EJA, and the real symmetric matrices. The others
517 will have to be overridden.
519 return (X
.conjugate_transpose()*Y
).trace()
525 Return the unit element of this algebra.
529 sage: from mjo.eja.eja_algebra import (HadamardEJA,
534 sage: J = HadamardEJA(5)
536 e0 + e1 + e2 + e3 + e4
540 The identity element acts like the identity::
542 sage: set_random_seed()
543 sage: J = random_eja()
544 sage: x = J.random_element()
545 sage: J.one()*x == x and x*J.one() == x
548 The matrix of the unit element's operator is the identity::
550 sage: set_random_seed()
551 sage: J = random_eja()
552 sage: actual = J.one().operator().matrix()
553 sage: expected = matrix.identity(J.base_ring(), J.dimension())
554 sage: actual == expected
558 # We can brute-force compute the matrices of the operators
559 # that correspond to the basis elements of this algebra.
560 # If some linear combination of those basis elements is the
561 # algebra identity, then the same linear combination of
562 # their matrices has to be the identity matrix.
564 # Of course, matrices aren't vectors in sage, so we have to
565 # appeal to the "long vectors" isometry.
566 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
568 # Now we use basis linear algebra to find the coefficients,
569 # of the matrices-as-vectors-linear-combination, which should
570 # work for the original algebra basis too.
571 A
= matrix
.column(self
.base_ring(), oper_vecs
)
573 # We used the isometry on the left-hand side already, but we
574 # still need to do it for the right-hand side. Recall that we
575 # wanted something that summed to the identity matrix.
576 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
578 # Now if there's an identity element in the algebra, this should work.
579 coeffs
= A
.solve_right(b
)
580 return self
.linear_combination(zip(self
.gens(), coeffs
))
583 def peirce_decomposition(self
, c
):
585 The Peirce decomposition of this algebra relative to the
588 In the future, this can be extended to a complete system of
589 orthogonal idempotents.
593 - ``c`` -- an idempotent of this algebra.
597 A triple (J0, J5, J1) containing two subalgebras and one subspace
600 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
601 corresponding to the eigenvalue zero.
603 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
604 corresponding to the eigenvalue one-half.
606 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
607 corresponding to the eigenvalue one.
609 These are the only possible eigenspaces for that operator, and this
610 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
611 orthogonal, and are subalgebras of this algebra with the appropriate
616 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
620 The canonical example comes from the symmetric matrices, which
621 decompose into diagonal and off-diagonal parts::
623 sage: J = RealSymmetricEJA(3)
624 sage: C = matrix(QQ, [ [1,0,0],
628 sage: J0,J5,J1 = J.peirce_decomposition(c)
630 Euclidean Jordan algebra of dimension 1...
632 Vector space of degree 6 and dimension 2...
634 Euclidean Jordan algebra of dimension 3...
638 Every algebra decomposes trivially with respect to its identity
641 sage: set_random_seed()
642 sage: J = random_eja()
643 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
644 sage: J0.dimension() == 0 and J5.dimension() == 0
646 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
649 The identity elements in the two subalgebras are the
650 projections onto their respective subspaces of the
651 superalgebra's identity element::
653 sage: set_random_seed()
654 sage: J = random_eja()
655 sage: x = J.random_element()
656 sage: if not J.is_trivial():
657 ....: while x.is_nilpotent():
658 ....: x = J.random_element()
659 sage: c = x.subalgebra_idempotent()
660 sage: J0,J5,J1 = J.peirce_decomposition(c)
661 sage: J1(c) == J1.one()
663 sage: J0(J.one() - c) == J0.one()
667 if not c
.is_idempotent():
668 raise ValueError("element is not idempotent: %s" % c
)
670 # Default these to what they should be if they turn out to be
671 # trivial, because eigenspaces_left() won't return eigenvalues
672 # corresponding to trivial spaces (e.g. it returns only the
673 # eigenspace corresponding to lambda=1 if you take the
674 # decomposition relative to the identity element).
675 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
676 J0
= trivial
# eigenvalue zero
677 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
678 J1
= trivial
# eigenvalue one
680 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
681 if eigval
== ~
(self
.base_ring()(2)):
684 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
685 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
691 raise ValueError("unexpected eigenvalue: %s" % eigval
)
696 def random_elements(self
, count
):
698 Return ``count`` random elements as a tuple.
702 sage: from mjo.eja.eja_algebra import JordanSpinEJA
706 sage: J = JordanSpinEJA(3)
707 sage: x,y,z = J.random_elements(3)
708 sage: all( [ x in J, y in J, z in J ])
710 sage: len( J.random_elements(10) ) == 10
714 return tuple( self
.random_element() for idx
in range(count
) )
717 def random_instance(cls
, field
=AA
, **kwargs
):
719 Return a random instance of this type of algebra.
721 Beware, this will crash for "most instances" because the
722 constructor below looks wrong.
724 if cls
is TrivialEJA
:
725 # The TrivialEJA class doesn't take an "n" argument because
729 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
730 return cls(n
, field
, **kwargs
)
733 def _charpoly_coefficients(self
):
735 The `r` polynomial coefficients of the "characteristic polynomial
739 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
740 R
= PolynomialRing(self
.base_ring(), var_names
)
742 F
= R
.fraction_field()
745 # From a result in my book, these are the entries of the
746 # basis representation of L_x.
747 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
750 L_x
= matrix(F
, n
, n
, L_x_i_j
)
751 # Compute an extra power in case the rank is equal to
752 # the dimension (otherwise, we would stop at x^(r-1)).
753 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
754 for k
in range(n
+1) ]
755 A
= matrix
.column(F
, x_powers
[:n
])
756 AE
= A
.extended_echelon_form()
762 # The theory says that only the first "r" coefficients are
763 # nonzero, and they actually live in the original polynomial
764 # ring and not the fraction field. We negate them because
765 # in the actual characteristic polynomial, they get moved
766 # to the other side where x^r lives.
767 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
772 Return the rank of this EJA.
774 This is a cached method because we know the rank a priori for
775 all of the algebras we can construct. Thus we can avoid the
776 expensive ``_charpoly_coefficients()`` call unless we truly
777 need to compute the whole characteristic polynomial.
781 sage: from mjo.eja.eja_algebra import (HadamardEJA,
783 ....: RealSymmetricEJA,
784 ....: ComplexHermitianEJA,
785 ....: QuaternionHermitianEJA,
790 The rank of the Jordan spin algebra is always two::
792 sage: JordanSpinEJA(2).rank()
794 sage: JordanSpinEJA(3).rank()
796 sage: JordanSpinEJA(4).rank()
799 The rank of the `n`-by-`n` Hermitian real, complex, or
800 quaternion matrices is `n`::
802 sage: RealSymmetricEJA(4).rank()
804 sage: ComplexHermitianEJA(3).rank()
806 sage: QuaternionHermitianEJA(2).rank()
811 Ensure that every EJA that we know how to construct has a
812 positive integer rank, unless the algebra is trivial in
813 which case its rank will be zero::
815 sage: set_random_seed()
816 sage: J = random_eja()
820 sage: r > 0 or (r == 0 and J.is_trivial())
823 Ensure that computing the rank actually works, since the ranks
824 of all simple algebras are known and will be cached by default::
826 sage: J = HadamardEJA(4)
827 sage: J.rank.clear_cache()
833 sage: J = JordanSpinEJA(4)
834 sage: J.rank.clear_cache()
840 sage: J = RealSymmetricEJA(3)
841 sage: J.rank.clear_cache()
847 sage: J = ComplexHermitianEJA(2)
848 sage: J.rank.clear_cache()
854 sage: J = QuaternionHermitianEJA(2)
855 sage: J.rank.clear_cache()
859 return len(self
._charpoly
_coefficients
())
862 def vector_space(self
):
864 Return the vector space that underlies this algebra.
868 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
872 sage: J = RealSymmetricEJA(2)
873 sage: J.vector_space()
874 Vector space of dimension 3 over...
877 return self
.zero().to_vector().parent().ambient_vector_space()
880 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
883 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
):
885 Return the Euclidean Jordan Algebra corresponding to the set
886 `R^n` under the Hadamard product.
888 Note: this is nothing more than the Cartesian product of ``n``
889 copies of the spin algebra. Once Cartesian product algebras
890 are implemented, this can go.
894 sage: from mjo.eja.eja_algebra import HadamardEJA
898 This multiplication table can be verified by hand::
900 sage: J = HadamardEJA(3)
901 sage: e0,e1,e2 = J.gens()
917 We can change the generator prefix::
919 sage: HadamardEJA(3, prefix='r').gens()
923 def __init__(self
, n
, field
=AA
, **kwargs
):
924 V
= VectorSpace(field
, n
)
925 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
928 fdeja
= super(HadamardEJA
, self
)
929 fdeja
.__init
__(field
, mult_table
, **kwargs
)
930 self
.rank
.set_cache(n
)
932 def inner_product(self
, x
, y
):
934 Faster to reimplement than to use natural representations.
938 sage: from mjo.eja.eja_algebra import HadamardEJA
942 Ensure that this is the usual inner product for the algebras
945 sage: set_random_seed()
946 sage: J = HadamardEJA.random_instance()
947 sage: x,y = J.random_elements(2)
948 sage: X = x.natural_representation()
949 sage: Y = y.natural_representation()
950 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
954 return x
.to_vector().inner_product(y
.to_vector())
957 def random_eja(field
=AA
, nontrivial
=False):
959 Return a "random" finite-dimensional Euclidean Jordan Algebra.
963 sage: from mjo.eja.eja_algebra import random_eja
968 Euclidean Jordan algebra of dimension...
971 eja_classes
= [HadamardEJA
,
975 QuaternionHermitianEJA
]
977 eja_classes
.append(TrivialEJA
)
978 classname
= choice(eja_classes
)
979 return classname
.random_instance(field
=field
)
986 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
988 def _max_test_case_size():
989 # Play it safe, since this will be squared and the underlying
990 # field can have dimension 4 (quaternions) too.
993 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
995 Compared to the superclass constructor, we take a basis instead of
996 a multiplication table because the latter can be computed in terms
997 of the former when the product is known (like it is here).
999 # Used in this class's fast _charpoly_coefficients() override.
1000 self
._basis
_normalizers
= None
1002 # We're going to loop through this a few times, so now's a good
1003 # time to ensure that it isn't a generator expression.
1004 basis
= tuple(basis
)
1006 if len(basis
) > 1 and normalize_basis
:
1007 # We'll need sqrt(2) to normalize the basis, and this
1008 # winds up in the multiplication table, so the whole
1009 # algebra needs to be over the field extension.
1010 R
= PolynomialRing(field
, 'z')
1013 if p
.is_irreducible():
1014 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1015 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1016 self
._basis
_normalizers
= tuple(
1017 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1018 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1020 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1022 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
1023 fdeja
.__init
__(field
, Qs
, natural_basis
=basis
, **kwargs
)
1028 def _charpoly_coefficients(self
):
1030 Override the parent method with something that tries to compute
1031 over a faster (non-extension) field.
1033 if self
._basis
_normalizers
is None:
1034 # We didn't normalize, so assume that the basis we started
1035 # with had entries in a nice field.
1036 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1038 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1039 self
._basis
_normalizers
) )
1041 # Do this over the rationals and convert back at the end.
1042 # Only works because we know the entries of the basis are
1044 J
= MatrixEuclideanJordanAlgebra(QQ
,
1046 normalize_basis
=False)
1047 return J
._charpoly
_coefficients
()
1051 def multiplication_table_from_matrix_basis(basis
):
1053 At least three of the five simple Euclidean Jordan algebras have the
1054 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1055 multiplication on the right is matrix multiplication. Given a basis
1056 for the underlying matrix space, this function returns a
1057 multiplication table (obtained by looping through the basis
1058 elements) for an algebra of those matrices.
1060 # In S^2, for example, we nominally have four coordinates even
1061 # though the space is of dimension three only. The vector space V
1062 # is supposed to hold the entire long vector, and the subspace W
1063 # of V will be spanned by the vectors that arise from symmetric
1064 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1065 field
= basis
[0].base_ring()
1066 dimension
= basis
[0].nrows()
1068 V
= VectorSpace(field
, dimension
**2)
1069 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1071 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1074 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1075 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1083 Embed the matrix ``M`` into a space of real matrices.
1085 The matrix ``M`` can have entries in any field at the moment:
1086 the real numbers, complex numbers, or quaternions. And although
1087 they are not a field, we can probably support octonions at some
1088 point, too. This function returns a real matrix that "acts like"
1089 the original with respect to matrix multiplication; i.e.
1091 real_embed(M*N) = real_embed(M)*real_embed(N)
1094 raise NotImplementedError
1098 def real_unembed(M
):
1100 The inverse of :meth:`real_embed`.
1102 raise NotImplementedError
1106 def natural_inner_product(cls
,X
,Y
):
1107 Xu
= cls
.real_unembed(X
)
1108 Yu
= cls
.real_unembed(Y
)
1109 tr
= (Xu
*Yu
).trace()
1112 # It's real already.
1115 # Otherwise, try the thing that works for complex numbers; and
1116 # if that doesn't work, the thing that works for quaternions.
1118 return tr
.vector()[0] # real part, imag part is index 1
1119 except AttributeError:
1120 # A quaternions doesn't have a vector() method, but does
1121 # have coefficient_tuple() method that returns the
1122 # coefficients of 1, i, j, and k -- in that order.
1123 return tr
.coefficient_tuple()[0]
1126 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1130 The identity function, for embedding real matrices into real
1136 def real_unembed(M
):
1138 The identity function, for unembedding real matrices from real
1144 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1146 The rank-n simple EJA consisting of real symmetric n-by-n
1147 matrices, the usual symmetric Jordan product, and the trace inner
1148 product. It has dimension `(n^2 + n)/2` over the reals.
1152 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1156 sage: J = RealSymmetricEJA(2)
1157 sage: e0, e1, e2 = J.gens()
1165 In theory, our "field" can be any subfield of the reals::
1167 sage: RealSymmetricEJA(2, RDF)
1168 Euclidean Jordan algebra of dimension 3 over Real Double Field
1169 sage: RealSymmetricEJA(2, RR)
1170 Euclidean Jordan algebra of dimension 3 over Real Field with
1171 53 bits of precision
1175 The dimension of this algebra is `(n^2 + n) / 2`::
1177 sage: set_random_seed()
1178 sage: n_max = RealSymmetricEJA._max_test_case_size()
1179 sage: n = ZZ.random_element(1, n_max)
1180 sage: J = RealSymmetricEJA(n)
1181 sage: J.dimension() == (n^2 + n)/2
1184 The Jordan multiplication is what we think it is::
1186 sage: set_random_seed()
1187 sage: J = RealSymmetricEJA.random_instance()
1188 sage: x,y = J.random_elements(2)
1189 sage: actual = (x*y).natural_representation()
1190 sage: X = x.natural_representation()
1191 sage: Y = y.natural_representation()
1192 sage: expected = (X*Y + Y*X)/2
1193 sage: actual == expected
1195 sage: J(expected) == x*y
1198 We can change the generator prefix::
1200 sage: RealSymmetricEJA(3, prefix='q').gens()
1201 (q0, q1, q2, q3, q4, q5)
1203 Our natural basis is normalized with respect to the natural inner
1204 product unless we specify otherwise::
1206 sage: set_random_seed()
1207 sage: J = RealSymmetricEJA.random_instance()
1208 sage: all( b.norm() == 1 for b in J.gens() )
1211 Since our natural basis is normalized with respect to the natural
1212 inner product, and since we know that this algebra is an EJA, any
1213 left-multiplication operator's matrix will be symmetric because
1214 natural->EJA basis representation is an isometry and within the EJA
1215 the operator is self-adjoint by the Jordan axiom::
1217 sage: set_random_seed()
1218 sage: x = RealSymmetricEJA.random_instance().random_element()
1219 sage: x.operator().matrix().is_symmetric()
1224 def _denormalized_basis(cls
, n
, field
):
1226 Return a basis for the space of real symmetric n-by-n matrices.
1230 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1234 sage: set_random_seed()
1235 sage: n = ZZ.random_element(1,5)
1236 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1237 sage: all( M.is_symmetric() for M in B)
1241 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1245 for j
in range(i
+1):
1246 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1250 Sij
= Eij
+ Eij
.transpose()
1256 def _max_test_case_size():
1257 return 4 # Dimension 10
1260 def __init__(self
, n
, field
=AA
, **kwargs
):
1261 basis
= self
._denormalized
_basis
(n
, field
)
1262 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, **kwargs
)
1263 self
.rank
.set_cache(n
)
1266 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1270 Embed the n-by-n complex matrix ``M`` into the space of real
1271 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1272 bi` to the block matrix ``[[a,b],[-b,a]]``.
1276 sage: from mjo.eja.eja_algebra import \
1277 ....: ComplexMatrixEuclideanJordanAlgebra
1281 sage: F = QuadraticField(-1, 'I')
1282 sage: x1 = F(4 - 2*i)
1283 sage: x2 = F(1 + 2*i)
1286 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1287 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1296 Embedding is a homomorphism (isomorphism, in fact)::
1298 sage: set_random_seed()
1299 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1300 sage: n = ZZ.random_element(n_max)
1301 sage: F = QuadraticField(-1, 'I')
1302 sage: X = random_matrix(F, n)
1303 sage: Y = random_matrix(F, n)
1304 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1305 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1306 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1313 raise ValueError("the matrix 'M' must be square")
1315 # We don't need any adjoined elements...
1316 field
= M
.base_ring().base_ring()
1320 a
= z
.list()[0] # real part, I guess
1321 b
= z
.list()[1] # imag part, I guess
1322 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1324 return matrix
.block(field
, n
, blocks
)
1328 def real_unembed(M
):
1330 The inverse of _embed_complex_matrix().
1334 sage: from mjo.eja.eja_algebra import \
1335 ....: ComplexMatrixEuclideanJordanAlgebra
1339 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1340 ....: [-2, 1, -4, 3],
1341 ....: [ 9, 10, 11, 12],
1342 ....: [-10, 9, -12, 11] ])
1343 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1345 [ 10*I + 9 12*I + 11]
1349 Unembedding is the inverse of embedding::
1351 sage: set_random_seed()
1352 sage: F = QuadraticField(-1, 'I')
1353 sage: M = random_matrix(F, 3)
1354 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1355 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1361 raise ValueError("the matrix 'M' must be square")
1362 if not n
.mod(2).is_zero():
1363 raise ValueError("the matrix 'M' must be a complex embedding")
1365 # If "M" was normalized, its base ring might have roots
1366 # adjoined and they can stick around after unembedding.
1367 field
= M
.base_ring()
1368 R
= PolynomialRing(field
, 'z')
1371 # Sage doesn't know how to embed AA into QQbar, i.e. how
1372 # to adjoin sqrt(-1) to AA.
1375 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1378 # Go top-left to bottom-right (reading order), converting every
1379 # 2-by-2 block we see to a single complex element.
1381 for k
in range(n
/2):
1382 for j
in range(n
/2):
1383 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1384 if submat
[0,0] != submat
[1,1]:
1385 raise ValueError('bad on-diagonal submatrix')
1386 if submat
[0,1] != -submat
[1,0]:
1387 raise ValueError('bad off-diagonal submatrix')
1388 z
= submat
[0,0] + submat
[0,1]*i
1391 return matrix(F
, n
/2, elements
)
1395 def natural_inner_product(cls
,X
,Y
):
1397 Compute a natural inner product in this algebra directly from
1402 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1406 This gives the same answer as the slow, default method implemented
1407 in :class:`MatrixEuclideanJordanAlgebra`::
1409 sage: set_random_seed()
1410 sage: J = ComplexHermitianEJA.random_instance()
1411 sage: x,y = J.random_elements(2)
1412 sage: Xe = x.natural_representation()
1413 sage: Ye = y.natural_representation()
1414 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1415 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1416 sage: expected = (X*Y).trace().real()
1417 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1418 sage: actual == expected
1422 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1425 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1427 The rank-n simple EJA consisting of complex Hermitian n-by-n
1428 matrices over the real numbers, the usual symmetric Jordan product,
1429 and the real-part-of-trace inner product. It has dimension `n^2` over
1434 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1438 In theory, our "field" can be any subfield of the reals::
1440 sage: ComplexHermitianEJA(2, RDF)
1441 Euclidean Jordan algebra of dimension 4 over Real Double Field
1442 sage: ComplexHermitianEJA(2, RR)
1443 Euclidean Jordan algebra of dimension 4 over Real Field with
1444 53 bits of precision
1448 The dimension of this algebra is `n^2`::
1450 sage: set_random_seed()
1451 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1452 sage: n = ZZ.random_element(1, n_max)
1453 sage: J = ComplexHermitianEJA(n)
1454 sage: J.dimension() == n^2
1457 The Jordan multiplication is what we think it is::
1459 sage: set_random_seed()
1460 sage: J = ComplexHermitianEJA.random_instance()
1461 sage: x,y = J.random_elements(2)
1462 sage: actual = (x*y).natural_representation()
1463 sage: X = x.natural_representation()
1464 sage: Y = y.natural_representation()
1465 sage: expected = (X*Y + Y*X)/2
1466 sage: actual == expected
1468 sage: J(expected) == x*y
1471 We can change the generator prefix::
1473 sage: ComplexHermitianEJA(2, prefix='z').gens()
1476 Our natural basis is normalized with respect to the natural inner
1477 product unless we specify otherwise::
1479 sage: set_random_seed()
1480 sage: J = ComplexHermitianEJA.random_instance()
1481 sage: all( b.norm() == 1 for b in J.gens() )
1484 Since our natural basis is normalized with respect to the natural
1485 inner product, and since we know that this algebra is an EJA, any
1486 left-multiplication operator's matrix will be symmetric because
1487 natural->EJA basis representation is an isometry and within the EJA
1488 the operator is self-adjoint by the Jordan axiom::
1490 sage: set_random_seed()
1491 sage: x = ComplexHermitianEJA.random_instance().random_element()
1492 sage: x.operator().matrix().is_symmetric()
1498 def _denormalized_basis(cls
, n
, field
):
1500 Returns a basis for the space of complex Hermitian n-by-n matrices.
1502 Why do we embed these? Basically, because all of numerical linear
1503 algebra assumes that you're working with vectors consisting of `n`
1504 entries from a field and scalars from the same field. There's no way
1505 to tell SageMath that (for example) the vectors contain complex
1506 numbers, while the scalar field is real.
1510 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1514 sage: set_random_seed()
1515 sage: n = ZZ.random_element(1,5)
1516 sage: field = QuadraticField(2, 'sqrt2')
1517 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1518 sage: all( M.is_symmetric() for M in B)
1522 R
= PolynomialRing(field
, 'z')
1524 F
= field
.extension(z
**2 + 1, 'I')
1527 # This is like the symmetric case, but we need to be careful:
1529 # * We want conjugate-symmetry, not just symmetry.
1530 # * The diagonal will (as a result) be real.
1534 for j
in range(i
+1):
1535 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1537 Sij
= cls
.real_embed(Eij
)
1540 # The second one has a minus because it's conjugated.
1541 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1543 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1546 # Since we embedded these, we can drop back to the "field" that we
1547 # started with instead of the complex extension "F".
1548 return ( s
.change_ring(field
) for s
in S
)
1551 def __init__(self
, n
, field
=AA
, **kwargs
):
1552 basis
= self
._denormalized
_basis
(n
,field
)
1553 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1554 self
.rank
.set_cache(n
)
1557 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1561 Embed the n-by-n quaternion matrix ``M`` into the space of real
1562 matrices of size 4n-by-4n by first sending each quaternion entry `z
1563 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1564 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1569 sage: from mjo.eja.eja_algebra import \
1570 ....: QuaternionMatrixEuclideanJordanAlgebra
1574 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1575 sage: i,j,k = Q.gens()
1576 sage: x = 1 + 2*i + 3*j + 4*k
1577 sage: M = matrix(Q, 1, [[x]])
1578 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1584 Embedding is a homomorphism (isomorphism, in fact)::
1586 sage: set_random_seed()
1587 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1588 sage: n = ZZ.random_element(n_max)
1589 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1590 sage: X = random_matrix(Q, n)
1591 sage: Y = random_matrix(Q, n)
1592 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1593 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1594 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1599 quaternions
= M
.base_ring()
1602 raise ValueError("the matrix 'M' must be square")
1604 F
= QuadraticField(-1, 'I')
1609 t
= z
.coefficient_tuple()
1614 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1615 [-c
+ d
*i
, a
- b
*i
]])
1616 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1617 blocks
.append(realM
)
1619 # We should have real entries by now, so use the realest field
1620 # we've got for the return value.
1621 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1626 def real_unembed(M
):
1628 The inverse of _embed_quaternion_matrix().
1632 sage: from mjo.eja.eja_algebra import \
1633 ....: QuaternionMatrixEuclideanJordanAlgebra
1637 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1638 ....: [-2, 1, -4, 3],
1639 ....: [-3, 4, 1, -2],
1640 ....: [-4, -3, 2, 1]])
1641 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1642 [1 + 2*i + 3*j + 4*k]
1646 Unembedding is the inverse of embedding::
1648 sage: set_random_seed()
1649 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1650 sage: M = random_matrix(Q, 3)
1651 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1652 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1658 raise ValueError("the matrix 'M' must be square")
1659 if not n
.mod(4).is_zero():
1660 raise ValueError("the matrix 'M' must be a quaternion embedding")
1662 # Use the base ring of the matrix to ensure that its entries can be
1663 # multiplied by elements of the quaternion algebra.
1664 field
= M
.base_ring()
1665 Q
= QuaternionAlgebra(field
,-1,-1)
1668 # Go top-left to bottom-right (reading order), converting every
1669 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1672 for l
in range(n
/4):
1673 for m
in range(n
/4):
1674 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1675 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1676 if submat
[0,0] != submat
[1,1].conjugate():
1677 raise ValueError('bad on-diagonal submatrix')
1678 if submat
[0,1] != -submat
[1,0].conjugate():
1679 raise ValueError('bad off-diagonal submatrix')
1680 z
= submat
[0,0].real()
1681 z
+= submat
[0,0].imag()*i
1682 z
+= submat
[0,1].real()*j
1683 z
+= submat
[0,1].imag()*k
1686 return matrix(Q
, n
/4, elements
)
1690 def natural_inner_product(cls
,X
,Y
):
1692 Compute a natural inner product in this algebra directly from
1697 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1701 This gives the same answer as the slow, default method implemented
1702 in :class:`MatrixEuclideanJordanAlgebra`::
1704 sage: set_random_seed()
1705 sage: J = QuaternionHermitianEJA.random_instance()
1706 sage: x,y = J.random_elements(2)
1707 sage: Xe = x.natural_representation()
1708 sage: Ye = y.natural_representation()
1709 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1710 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1711 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1712 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1713 sage: actual == expected
1717 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1720 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1722 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1723 matrices, the usual symmetric Jordan product, and the
1724 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1729 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1733 In theory, our "field" can be any subfield of the reals::
1735 sage: QuaternionHermitianEJA(2, RDF)
1736 Euclidean Jordan algebra of dimension 6 over Real Double Field
1737 sage: QuaternionHermitianEJA(2, RR)
1738 Euclidean Jordan algebra of dimension 6 over Real Field with
1739 53 bits of precision
1743 The dimension of this algebra is `2*n^2 - n`::
1745 sage: set_random_seed()
1746 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1747 sage: n = ZZ.random_element(1, n_max)
1748 sage: J = QuaternionHermitianEJA(n)
1749 sage: J.dimension() == 2*(n^2) - n
1752 The Jordan multiplication is what we think it is::
1754 sage: set_random_seed()
1755 sage: J = QuaternionHermitianEJA.random_instance()
1756 sage: x,y = J.random_elements(2)
1757 sage: actual = (x*y).natural_representation()
1758 sage: X = x.natural_representation()
1759 sage: Y = y.natural_representation()
1760 sage: expected = (X*Y + Y*X)/2
1761 sage: actual == expected
1763 sage: J(expected) == x*y
1766 We can change the generator prefix::
1768 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1769 (a0, a1, a2, a3, a4, a5)
1771 Our natural basis is normalized with respect to the natural inner
1772 product unless we specify otherwise::
1774 sage: set_random_seed()
1775 sage: J = QuaternionHermitianEJA.random_instance()
1776 sage: all( b.norm() == 1 for b in J.gens() )
1779 Since our natural basis is normalized with respect to the natural
1780 inner product, and since we know that this algebra is an EJA, any
1781 left-multiplication operator's matrix will be symmetric because
1782 natural->EJA basis representation is an isometry and within the EJA
1783 the operator is self-adjoint by the Jordan axiom::
1785 sage: set_random_seed()
1786 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1787 sage: x.operator().matrix().is_symmetric()
1792 def _denormalized_basis(cls
, n
, field
):
1794 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1796 Why do we embed these? Basically, because all of numerical
1797 linear algebra assumes that you're working with vectors consisting
1798 of `n` entries from a field and scalars from the same field. There's
1799 no way to tell SageMath that (for example) the vectors contain
1800 complex numbers, while the scalar field is real.
1804 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1808 sage: set_random_seed()
1809 sage: n = ZZ.random_element(1,5)
1810 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1811 sage: all( M.is_symmetric() for M in B )
1815 Q
= QuaternionAlgebra(QQ
,-1,-1)
1818 # This is like the symmetric case, but we need to be careful:
1820 # * We want conjugate-symmetry, not just symmetry.
1821 # * The diagonal will (as a result) be real.
1825 for j
in range(i
+1):
1826 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1828 Sij
= cls
.real_embed(Eij
)
1831 # The second, third, and fourth ones have a minus
1832 # because they're conjugated.
1833 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1835 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1837 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1839 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1842 # Since we embedded these, we can drop back to the "field" that we
1843 # started with instead of the quaternion algebra "Q".
1844 return ( s
.change_ring(field
) for s
in S
)
1847 def __init__(self
, n
, field
=AA
, **kwargs
):
1848 basis
= self
._denormalized
_basis
(n
,field
)
1849 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1850 self
.rank
.set_cache(n
)
1853 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1855 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1856 with the half-trace inner product and jordan product ``x*y =
1857 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
1858 symmetric positive-definite "bilinear form" matrix. It has
1859 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
1860 when ``B`` is the identity matrix of order ``n-1``.
1864 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1865 ....: JordanSpinEJA)
1869 When no bilinear form is specified, the identity matrix is used,
1870 and the resulting algebra is the Jordan spin algebra::
1872 sage: J0 = BilinearFormEJA(3)
1873 sage: J1 = JordanSpinEJA(3)
1874 sage: J0.multiplication_table() == J0.multiplication_table()
1879 We can create a zero-dimensional algebra::
1881 sage: J = BilinearFormEJA(0)
1885 We can check the multiplication condition given in the Jordan, von
1886 Neumann, and Wigner paper (and also discussed on my "On the
1887 symmetry..." paper). Note that this relies heavily on the standard
1888 choice of basis, as does anything utilizing the bilinear form matrix::
1890 sage: set_random_seed()
1891 sage: n = ZZ.random_element(5)
1892 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
1893 sage: B = M.transpose()*M
1894 sage: J = BilinearFormEJA(n, B=B)
1895 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
1896 sage: V = J.vector_space()
1897 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
1898 ....: for ei in eis ]
1899 sage: actual = [ sis[i]*sis[j]
1900 ....: for i in range(n-1)
1901 ....: for j in range(n-1) ]
1902 sage: expected = [ J.one() if i == j else J.zero()
1903 ....: for i in range(n-1)
1904 ....: for j in range(n-1) ]
1905 sage: actual == expected
1908 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
1910 self
._B
= matrix
.identity(field
, max(0,n
-1))
1914 V
= VectorSpace(field
, n
)
1915 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1924 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
1925 zbar
= y0
*xbar
+ x0
*ybar
1926 z
= V([z0
] + zbar
.list())
1927 mult_table
[i
][j
] = z
1929 # The rank of this algebra is two, unless we're in a
1930 # one-dimensional ambient space (because the rank is bounded
1931 # by the ambient dimension).
1932 fdeja
= super(BilinearFormEJA
, self
)
1933 fdeja
.__init
__(field
, mult_table
, **kwargs
)
1934 self
.rank
.set_cache(min(n
,2))
1936 def inner_product(self
, x
, y
):
1938 Half of the trace inner product.
1940 This is defined so that the special case of the Jordan spin
1941 algebra gets the usual inner product.
1945 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1949 Ensure that this is one-half of the trace inner-product when
1950 the algebra isn't just the reals (when ``n`` isn't one). This
1951 is in Faraut and Koranyi, and also my "On the symmetry..."
1954 sage: set_random_seed()
1955 sage: n = ZZ.random_element(2,5)
1956 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
1957 sage: B = M.transpose()*M
1958 sage: J = BilinearFormEJA(n, B=B)
1959 sage: x = J.random_element()
1960 sage: y = J.random_element()
1961 sage: x.inner_product(y) == (x*y).trace()/2
1965 xvec
= x
.to_vector()
1967 yvec
= y
.to_vector()
1969 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
1972 class JordanSpinEJA(BilinearFormEJA
):
1974 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1975 with the usual inner product and jordan product ``x*y =
1976 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1981 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1985 This multiplication table can be verified by hand::
1987 sage: J = JordanSpinEJA(4)
1988 sage: e0,e1,e2,e3 = J.gens()
2004 We can change the generator prefix::
2006 sage: JordanSpinEJA(2, prefix='B').gens()
2011 Ensure that we have the usual inner product on `R^n`::
2013 sage: set_random_seed()
2014 sage: J = JordanSpinEJA.random_instance()
2015 sage: x,y = J.random_elements(2)
2016 sage: X = x.natural_representation()
2017 sage: Y = y.natural_representation()
2018 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2022 def __init__(self
, n
, field
=AA
, **kwargs
):
2023 # This is a special case of the BilinearFormEJA with the identity
2024 # matrix as its bilinear form.
2025 return super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2028 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2030 The trivial Euclidean Jordan algebra consisting of only a zero element.
2034 sage: from mjo.eja.eja_algebra import TrivialEJA
2038 sage: J = TrivialEJA()
2045 sage: 7*J.one()*12*J.one()
2047 sage: J.one().inner_product(J.one())
2049 sage: J.one().norm()
2051 sage: J.one().subalgebra_generated_by()
2052 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2057 def __init__(self
, field
=AA
, **kwargs
):
2059 fdeja
= super(TrivialEJA
, self
)
2060 # The rank is zero using my definition, namely the dimension of the
2061 # largest subalgebra generated by any element.
2062 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2063 self
.rank
.set_cache(0)