2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 from itertools
import repeat
10 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
11 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
12 from sage
.combinat
.free_module
import CombinatorialFreeModule
13 from sage
.matrix
.constructor
import matrix
14 from sage
.matrix
.matrix_space
import MatrixSpace
15 from sage
.misc
.cachefunc
import cached_method
16 from sage
.misc
.lazy_import
import lazy_import
17 from sage
.misc
.prandom
import choice
18 from sage
.misc
.table
import table
19 from sage
.modules
.free_module
import FreeModule
, VectorSpace
20 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
23 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
24 lazy_import('mjo.eja.eja_subalgebra',
25 'FiniteDimensionalEuclideanJordanSubalgebra')
26 from mjo
.eja
.eja_utils
import _mat2vec
28 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
30 def _coerce_map_from_base_ring(self
):
32 Disable the map from the base ring into the algebra.
34 Performing a nonsense conversion like this automatically
35 is counterpedagogical. The fallback is to try the usual
36 element constructor, which should also fail.
40 sage: from mjo.eja.eja_algebra import random_eja
44 sage: set_random_seed()
45 sage: J = random_eja()
47 Traceback (most recent call last):
49 ValueError: not a naturally-represented algebra element
64 sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
68 By definition, Jordan multiplication commutes::
70 sage: set_random_seed()
71 sage: J = random_eja()
72 sage: x,y = J.random_elements(2)
78 The ``field`` we're given must be real::
80 sage: JordanSpinEJA(2,QQbar)
81 Traceback (most recent call last):
83 ValueError: field is not real
87 if not field
.is_subring(RR
):
88 # Note: this does return true for the real algebraic
89 # field, and any quadratic field where we've specified
91 raise ValueError('field is not real')
93 self
._natural
_basis
= natural_basis
96 category
= MagmaticAlgebras(field
).FiniteDimensional()
97 category
= category
.WithBasis().Unital()
99 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
101 range(len(mult_table
)),
104 self
.print_options(bracket
='')
106 # The multiplication table we're given is necessarily in terms
107 # of vectors, because we don't have an algebra yet for
108 # anything to be an element of. However, it's faster in the
109 # long run to have the multiplication table be in terms of
110 # algebra elements. We do this after calling the superclass
111 # constructor so that from_vector() knows what to do.
112 self
._multiplication
_table
= [
113 list(map(lambda x
: self
.from_vector(x
), ls
))
118 def _element_constructor_(self
, elt
):
120 Construct an element of this algebra from its natural
123 This gets called only after the parent element _call_ method
124 fails to find a coercion for the argument.
128 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
130 ....: RealSymmetricEJA)
134 The identity in `S^n` is converted to the identity in the EJA::
136 sage: J = RealSymmetricEJA(3)
137 sage: I = matrix.identity(QQ,3)
138 sage: J(I) == J.one()
141 This skew-symmetric matrix can't be represented in the EJA::
143 sage: J = RealSymmetricEJA(3)
144 sage: A = matrix(QQ,3, lambda i,j: i-j)
146 Traceback (most recent call last):
148 ArithmeticError: vector is not in free module
152 Ensure that we can convert any element of the two non-matrix
153 simple algebras (whose natural representations are their usual
154 vector representations) back and forth faithfully::
156 sage: set_random_seed()
157 sage: J = HadamardEJA.random_instance()
158 sage: x = J.random_element()
159 sage: J(x.to_vector().column()) == x
161 sage: J = JordanSpinEJA.random_instance()
162 sage: x = J.random_element()
163 sage: J(x.to_vector().column()) == x
167 msg
= "not a naturally-represented algebra element"
169 # The superclass implementation of random_element()
170 # needs to be able to coerce "0" into the algebra.
172 elif elt
in self
.base_ring():
173 # Ensure that no base ring -> algebra coercion is performed
174 # by this method. There's some stupidity in sage that would
175 # otherwise propagate to this method; for example, sage thinks
176 # that the integer 3 belongs to the space of 2-by-2 matrices.
177 raise ValueError(msg
)
179 natural_basis
= self
.natural_basis()
180 basis_space
= natural_basis
[0].matrix_space()
181 if elt
not in basis_space
:
182 raise ValueError(msg
)
184 # Thanks for nothing! Matrix spaces aren't vector spaces in
185 # Sage, so we have to figure out its natural-basis coordinates
186 # ourselves. We use the basis space's ring instead of the
187 # element's ring because the basis space might be an algebraic
188 # closure whereas the base ring of the 3-by-3 identity matrix
189 # could be QQ instead of QQbar.
190 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
191 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
192 coords
= W
.coordinate_vector(_mat2vec(elt
))
193 return self
.from_vector(coords
)
196 def _max_test_case_size():
198 Return an integer "size" that is an upper bound on the size of
199 this algebra when it is used in a random test
200 case. Unfortunately, the term "size" is quite vague -- when
201 dealing with `R^n` under either the Hadamard or Jordan spin
202 product, the "size" refers to the dimension `n`. When dealing
203 with a matrix algebra (real symmetric or complex/quaternion
204 Hermitian), it refers to the size of the matrix, which is
205 far less than the dimension of the underlying vector space.
207 We default to five in this class, which is safe in `R^n`. The
208 matrix algebra subclasses (or any class where the "size" is
209 interpreted to be far less than the dimension) should override
210 with a smaller number.
216 Return a string representation of ``self``.
220 sage: from mjo.eja.eja_algebra import JordanSpinEJA
224 Ensure that it says what we think it says::
226 sage: JordanSpinEJA(2, field=AA)
227 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
228 sage: JordanSpinEJA(3, field=RDF)
229 Euclidean Jordan algebra of dimension 3 over Real Double Field
232 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
233 return fmt
.format(self
.dimension(), self
.base_ring())
235 def product_on_basis(self
, i
, j
):
236 return self
._multiplication
_table
[i
][j
]
239 def characteristic_polynomial(self
):
241 Return a characteristic polynomial that works for all elements
244 The resulting polynomial has `n+1` variables, where `n` is the
245 dimension of this algebra. The first `n` variables correspond to
246 the coordinates of an algebra element: when evaluated at the
247 coordinates of an algebra element with respect to a certain
248 basis, the result is a univariate polynomial (in the one
249 remaining variable ``t``), namely the characteristic polynomial
254 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
258 The characteristic polynomial in the spin algebra is given in
259 Alizadeh, Example 11.11::
261 sage: J = JordanSpinEJA(3)
262 sage: p = J.characteristic_polynomial(); p
263 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
264 sage: xvec = J.one().to_vector()
268 By definition, the characteristic polynomial is a monic
269 degree-zero polynomial in a rank-zero algebra. Note that
270 Cayley-Hamilton is indeed satisfied since the polynomial
271 ``1`` evaluates to the identity element of the algebra on
274 sage: J = TrivialEJA()
275 sage: J.characteristic_polynomial()
282 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
283 a
= self
._charpoly
_coefficients
()
285 # We go to a bit of trouble here to reorder the
286 # indeterminates, so that it's easier to evaluate the
287 # characteristic polynomial at x's coordinates and get back
288 # something in terms of t, which is what we want.
289 S
= PolynomialRing(self
.base_ring(),'t')
293 S
= PolynomialRing(S
, R
.variable_names())
296 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
299 def inner_product(self
, x
, y
):
301 The inner product associated with this Euclidean Jordan algebra.
303 Defaults to the trace inner product, but can be overridden by
304 subclasses if they are sure that the necessary properties are
309 sage: from mjo.eja.eja_algebra import random_eja
313 Our inner product is "associative," which means the following for
314 a symmetric bilinear form::
316 sage: set_random_seed()
317 sage: J = random_eja()
318 sage: x,y,z = J.random_elements(3)
319 sage: (x*y).inner_product(z) == y.inner_product(x*z)
323 X
= x
.natural_representation()
324 Y
= y
.natural_representation()
325 return self
.natural_inner_product(X
,Y
)
328 def is_trivial(self
):
330 Return whether or not this algebra is trivial.
332 A trivial algebra contains only the zero element.
336 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
341 sage: J = ComplexHermitianEJA(3)
347 sage: J = TrivialEJA()
352 return self
.dimension() == 0
355 def multiplication_table(self
):
357 Return a visual representation of this algebra's multiplication
358 table (on basis elements).
362 sage: from mjo.eja.eja_algebra import JordanSpinEJA
366 sage: J = JordanSpinEJA(4)
367 sage: J.multiplication_table()
368 +----++----+----+----+----+
369 | * || e0 | e1 | e2 | e3 |
370 +====++====+====+====+====+
371 | e0 || e0 | e1 | e2 | e3 |
372 +----++----+----+----+----+
373 | e1 || e1 | e0 | 0 | 0 |
374 +----++----+----+----+----+
375 | e2 || e2 | 0 | e0 | 0 |
376 +----++----+----+----+----+
377 | e3 || e3 | 0 | 0 | e0 |
378 +----++----+----+----+----+
381 M
= list(self
._multiplication
_table
) # copy
382 for i
in range(len(M
)):
383 # M had better be "square"
384 M
[i
] = [self
.monomial(i
)] + M
[i
]
385 M
= [["*"] + list(self
.gens())] + M
386 return table(M
, header_row
=True, header_column
=True, frame
=True)
389 def natural_basis(self
):
391 Return a more-natural representation of this algebra's basis.
393 Every finite-dimensional Euclidean Jordan Algebra is a direct
394 sum of five simple algebras, four of which comprise Hermitian
395 matrices. This method returns the original "natural" basis
396 for our underlying vector space. (Typically, the natural basis
397 is used to construct the multiplication table in the first place.)
399 Note that this will always return a matrix. The standard basis
400 in `R^n` will be returned as `n`-by-`1` column matrices.
404 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
405 ....: RealSymmetricEJA)
409 sage: J = RealSymmetricEJA(2)
411 Finite family {0: e0, 1: e1, 2: e2}
412 sage: J.natural_basis()
414 [1 0] [ 0 0.7071067811865475?] [0 0]
415 [0 0], [0.7071067811865475? 0], [0 1]
420 sage: J = JordanSpinEJA(2)
422 Finite family {0: e0, 1: e1}
423 sage: J.natural_basis()
430 if self
._natural
_basis
is None:
431 M
= self
.natural_basis_space()
432 return tuple( M(b
.to_vector()) for b
in self
.basis() )
434 return self
._natural
_basis
437 def natural_basis_space(self
):
439 Return the matrix space in which this algebra's natural basis
442 Generally this will be an `n`-by-`1` column-vector space,
443 except when the algebra is trivial. There it's `n`-by-`n`
444 (where `n` is zero), to ensure that two elements of the
445 natural basis space (empty matrices) can be multiplied.
447 if self
.is_trivial():
448 return MatrixSpace(self
.base_ring(), 0)
449 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
450 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
452 return self
._natural
_basis
[0].matrix_space()
456 def natural_inner_product(X
,Y
):
458 Compute the inner product of two naturally-represented elements.
460 For example in the real symmetric matrix EJA, this will compute
461 the trace inner-product of two n-by-n symmetric matrices. The
462 default should work for the real cartesian product EJA, the
463 Jordan spin EJA, and the real symmetric matrices. The others
464 will have to be overridden.
466 return (X
.conjugate_transpose()*Y
).trace()
472 Return the unit element of this algebra.
476 sage: from mjo.eja.eja_algebra import (HadamardEJA,
481 sage: J = HadamardEJA(5)
483 e0 + e1 + e2 + e3 + e4
487 The identity element acts like the identity::
489 sage: set_random_seed()
490 sage: J = random_eja()
491 sage: x = J.random_element()
492 sage: J.one()*x == x and x*J.one() == x
495 The matrix of the unit element's operator is the identity::
497 sage: set_random_seed()
498 sage: J = random_eja()
499 sage: actual = J.one().operator().matrix()
500 sage: expected = matrix.identity(J.base_ring(), J.dimension())
501 sage: actual == expected
505 # We can brute-force compute the matrices of the operators
506 # that correspond to the basis elements of this algebra.
507 # If some linear combination of those basis elements is the
508 # algebra identity, then the same linear combination of
509 # their matrices has to be the identity matrix.
511 # Of course, matrices aren't vectors in sage, so we have to
512 # appeal to the "long vectors" isometry.
513 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
515 # Now we use basis linear algebra to find the coefficients,
516 # of the matrices-as-vectors-linear-combination, which should
517 # work for the original algebra basis too.
518 A
= matrix
.column(self
.base_ring(), oper_vecs
)
520 # We used the isometry on the left-hand side already, but we
521 # still need to do it for the right-hand side. Recall that we
522 # wanted something that summed to the identity matrix.
523 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
525 # Now if there's an identity element in the algebra, this should work.
526 coeffs
= A
.solve_right(b
)
527 return self
.linear_combination(zip(self
.gens(), coeffs
))
530 def peirce_decomposition(self
, c
):
532 The Peirce decomposition of this algebra relative to the
535 In the future, this can be extended to a complete system of
536 orthogonal idempotents.
540 - ``c`` -- an idempotent of this algebra.
544 A triple (J0, J5, J1) containing two subalgebras and one subspace
547 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
548 corresponding to the eigenvalue zero.
550 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
551 corresponding to the eigenvalue one-half.
553 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
554 corresponding to the eigenvalue one.
556 These are the only possible eigenspaces for that operator, and this
557 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
558 orthogonal, and are subalgebras of this algebra with the appropriate
563 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
567 The canonical example comes from the symmetric matrices, which
568 decompose into diagonal and off-diagonal parts::
570 sage: J = RealSymmetricEJA(3)
571 sage: C = matrix(QQ, [ [1,0,0],
575 sage: J0,J5,J1 = J.peirce_decomposition(c)
577 Euclidean Jordan algebra of dimension 1...
579 Vector space of degree 6 and dimension 2...
581 Euclidean Jordan algebra of dimension 3...
585 Every algebra decomposes trivially with respect to its identity
588 sage: set_random_seed()
589 sage: J = random_eja()
590 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
591 sage: J0.dimension() == 0 and J5.dimension() == 0
593 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
596 The identity elements in the two subalgebras are the
597 projections onto their respective subspaces of the
598 superalgebra's identity element::
600 sage: set_random_seed()
601 sage: J = random_eja()
602 sage: x = J.random_element()
603 sage: if not J.is_trivial():
604 ....: while x.is_nilpotent():
605 ....: x = J.random_element()
606 sage: c = x.subalgebra_idempotent()
607 sage: J0,J5,J1 = J.peirce_decomposition(c)
608 sage: J1(c) == J1.one()
610 sage: J0(J.one() - c) == J0.one()
614 if not c
.is_idempotent():
615 raise ValueError("element is not idempotent: %s" % c
)
617 # Default these to what they should be if they turn out to be
618 # trivial, because eigenspaces_left() won't return eigenvalues
619 # corresponding to trivial spaces (e.g. it returns only the
620 # eigenspace corresponding to lambda=1 if you take the
621 # decomposition relative to the identity element).
622 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
623 J0
= trivial
# eigenvalue zero
624 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
625 J1
= trivial
# eigenvalue one
627 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
628 if eigval
== ~
(self
.base_ring()(2)):
631 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
632 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
638 raise ValueError("unexpected eigenvalue: %s" % eigval
)
643 def random_elements(self
, count
):
645 Return ``count`` random elements as a tuple.
649 sage: from mjo.eja.eja_algebra import JordanSpinEJA
653 sage: J = JordanSpinEJA(3)
654 sage: x,y,z = J.random_elements(3)
655 sage: all( [ x in J, y in J, z in J ])
657 sage: len( J.random_elements(10) ) == 10
661 return tuple( self
.random_element() for idx
in range(count
) )
664 def random_instance(cls
, field
=AA
, **kwargs
):
666 Return a random instance of this type of algebra.
668 Beware, this will crash for "most instances" because the
669 constructor below looks wrong.
671 if cls
is TrivialEJA
:
672 # The TrivialEJA class doesn't take an "n" argument because
676 n
= ZZ
.random_element(cls
._max
_test
_case
_size
() + 1)
677 return cls(n
, field
, **kwargs
)
680 def _charpoly_coefficients(self
):
682 The `r` polynomial coefficients of the "characteristic polynomial
686 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
687 R
= PolynomialRing(self
.base_ring(), var_names
)
689 F
= R
.fraction_field()
692 # From a result in my book, these are the entries of the
693 # basis representation of L_x.
694 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
697 L_x
= matrix(F
, n
, n
, L_x_i_j
)
700 if self
.rank
.is_in_cache():
702 # There's no need to pad the system with redundant
703 # columns if we *know* they'll be redundant.
706 # Compute an extra power in case the rank is equal to
707 # the dimension (otherwise, we would stop at x^(r-1)).
708 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
709 for k
in range(n
+1) ]
710 A
= matrix
.column(F
, x_powers
[:n
])
711 AE
= A
.extended_echelon_form()
718 # The theory says that only the first "r" coefficients are
719 # nonzero, and they actually live in the original polynomial
720 # ring and not the fraction field. We negate them because
721 # in the actual characteristic polynomial, they get moved
722 # to the other side where x^r lives.
723 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
728 Return the rank of this EJA.
730 This is a cached method because we know the rank a priori for
731 all of the algebras we can construct. Thus we can avoid the
732 expensive ``_charpoly_coefficients()`` call unless we truly
733 need to compute the whole characteristic polynomial.
737 sage: from mjo.eja.eja_algebra import (HadamardEJA,
739 ....: RealSymmetricEJA,
740 ....: ComplexHermitianEJA,
741 ....: QuaternionHermitianEJA,
746 The rank of the Jordan spin algebra is always two::
748 sage: JordanSpinEJA(2).rank()
750 sage: JordanSpinEJA(3).rank()
752 sage: JordanSpinEJA(4).rank()
755 The rank of the `n`-by-`n` Hermitian real, complex, or
756 quaternion matrices is `n`::
758 sage: RealSymmetricEJA(4).rank()
760 sage: ComplexHermitianEJA(3).rank()
762 sage: QuaternionHermitianEJA(2).rank()
767 Ensure that every EJA that we know how to construct has a
768 positive integer rank, unless the algebra is trivial in
769 which case its rank will be zero::
771 sage: set_random_seed()
772 sage: J = random_eja()
776 sage: r > 0 or (r == 0 and J.is_trivial())
779 Ensure that computing the rank actually works, since the ranks
780 of all simple algebras are known and will be cached by default::
782 sage: J = HadamardEJA(4)
783 sage: J.rank.clear_cache()
789 sage: J = JordanSpinEJA(4)
790 sage: J.rank.clear_cache()
796 sage: J = RealSymmetricEJA(3)
797 sage: J.rank.clear_cache()
803 sage: J = ComplexHermitianEJA(2)
804 sage: J.rank.clear_cache()
810 sage: J = QuaternionHermitianEJA(2)
811 sage: J.rank.clear_cache()
815 return len(self
._charpoly
_coefficients
())
818 def vector_space(self
):
820 Return the vector space that underlies this algebra.
824 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
828 sage: J = RealSymmetricEJA(2)
829 sage: J.vector_space()
830 Vector space of dimension 3 over...
833 return self
.zero().to_vector().parent().ambient_vector_space()
836 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
839 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
):
841 Return the Euclidean Jordan Algebra corresponding to the set
842 `R^n` under the Hadamard product.
844 Note: this is nothing more than the Cartesian product of ``n``
845 copies of the spin algebra. Once Cartesian product algebras
846 are implemented, this can go.
850 sage: from mjo.eja.eja_algebra import HadamardEJA
854 This multiplication table can be verified by hand::
856 sage: J = HadamardEJA(3)
857 sage: e0,e1,e2 = J.gens()
873 We can change the generator prefix::
875 sage: HadamardEJA(3, prefix='r').gens()
879 def __init__(self
, n
, field
=AA
, **kwargs
):
880 V
= VectorSpace(field
, n
)
881 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
884 fdeja
= super(HadamardEJA
, self
)
885 fdeja
.__init
__(field
, mult_table
, **kwargs
)
886 self
.rank
.set_cache(n
)
888 def inner_product(self
, x
, y
):
890 Faster to reimplement than to use natural representations.
894 sage: from mjo.eja.eja_algebra import HadamardEJA
898 Ensure that this is the usual inner product for the algebras
901 sage: set_random_seed()
902 sage: J = HadamardEJA.random_instance()
903 sage: x,y = J.random_elements(2)
904 sage: X = x.natural_representation()
905 sage: Y = y.natural_representation()
906 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
910 return x
.to_vector().inner_product(y
.to_vector())
913 def random_eja(field
=AA
):
915 Return a "random" finite-dimensional Euclidean Jordan Algebra.
919 sage: from mjo.eja.eja_algebra import random_eja
924 Euclidean Jordan algebra of dimension...
927 classname
= choice([TrivialEJA
,
932 QuaternionHermitianEJA
])
933 return classname
.random_instance(field
=field
)
938 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
940 def _max_test_case_size():
941 # Play it safe, since this will be squared and the underlying
942 # field can have dimension 4 (quaternions) too.
945 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
947 Compared to the superclass constructor, we take a basis instead of
948 a multiplication table because the latter can be computed in terms
949 of the former when the product is known (like it is here).
951 # Used in this class's fast _charpoly_coefficients() override.
952 self
._basis
_normalizers
= None
954 # We're going to loop through this a few times, so now's a good
955 # time to ensure that it isn't a generator expression.
958 if len(basis
) > 1 and normalize_basis
:
959 # We'll need sqrt(2) to normalize the basis, and this
960 # winds up in the multiplication table, so the whole
961 # algebra needs to be over the field extension.
962 R
= PolynomialRing(field
, 'z')
965 if p
.is_irreducible():
966 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
967 basis
= tuple( s
.change_ring(field
) for s
in basis
)
968 self
._basis
_normalizers
= tuple(
969 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
970 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
972 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
974 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
975 fdeja
.__init
__(field
, Qs
, natural_basis
=basis
, **kwargs
)
980 def _charpoly_coefficients(self
):
982 Override the parent method with something that tries to compute
983 over a faster (non-extension) field.
985 if self
._basis
_normalizers
is None:
986 # We didn't normalize, so assume that the basis we started
987 # with had entries in a nice field.
988 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
990 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
991 self
._basis
_normalizers
) )
993 # Do this over the rationals and convert back at the end.
994 # Only works because we know the entries of the basis are
996 J
= MatrixEuclideanJordanAlgebra(QQ
,
998 normalize_basis
=False)
999 a
= J
._charpoly
_coefficients
()
1001 # Unfortunately, changing the basis does change the
1002 # coefficients of the characteristic polynomial, but since
1003 # these are really the coefficients of the "characteristic
1004 # polynomial of" function, everything is still nice and
1005 # unevaluated. It's therefore "obvious" how scaling the
1006 # basis affects the coordinate variables X1, X2, et
1007 # cetera. Scaling the first basis vector up by "n" adds a
1008 # factor of 1/n into every "X1" term, for example. So here
1009 # we simply undo the basis_normalizer scaling that we
1010 # performed earlier.
1012 # The a[0] access here is safe because trivial algebras
1013 # won't have any basis normalizers and therefore won't
1014 # make it to this "else" branch.
1015 XS
= a
[0].parent().gens()
1016 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1017 for i
in range(len(XS
)) }
1018 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1022 def multiplication_table_from_matrix_basis(basis
):
1024 At least three of the five simple Euclidean Jordan algebras have the
1025 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1026 multiplication on the right is matrix multiplication. Given a basis
1027 for the underlying matrix space, this function returns a
1028 multiplication table (obtained by looping through the basis
1029 elements) for an algebra of those matrices.
1031 # In S^2, for example, we nominally have four coordinates even
1032 # though the space is of dimension three only. The vector space V
1033 # is supposed to hold the entire long vector, and the subspace W
1034 # of V will be spanned by the vectors that arise from symmetric
1035 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1039 field
= basis
[0].base_ring()
1040 dimension
= basis
[0].nrows()
1042 V
= VectorSpace(field
, dimension
**2)
1043 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1045 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1048 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1049 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1057 Embed the matrix ``M`` into a space of real matrices.
1059 The matrix ``M`` can have entries in any field at the moment:
1060 the real numbers, complex numbers, or quaternions. And although
1061 they are not a field, we can probably support octonions at some
1062 point, too. This function returns a real matrix that "acts like"
1063 the original with respect to matrix multiplication; i.e.
1065 real_embed(M*N) = real_embed(M)*real_embed(N)
1068 raise NotImplementedError
1072 def real_unembed(M
):
1074 The inverse of :meth:`real_embed`.
1076 raise NotImplementedError
1080 def natural_inner_product(cls
,X
,Y
):
1081 Xu
= cls
.real_unembed(X
)
1082 Yu
= cls
.real_unembed(Y
)
1083 tr
= (Xu
*Yu
).trace()
1086 # It's real already.
1089 # Otherwise, try the thing that works for complex numbers; and
1090 # if that doesn't work, the thing that works for quaternions.
1092 return tr
.vector()[0] # real part, imag part is index 1
1093 except AttributeError:
1094 # A quaternions doesn't have a vector() method, but does
1095 # have coefficient_tuple() method that returns the
1096 # coefficients of 1, i, j, and k -- in that order.
1097 return tr
.coefficient_tuple()[0]
1100 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1104 The identity function, for embedding real matrices into real
1110 def real_unembed(M
):
1112 The identity function, for unembedding real matrices from real
1118 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1120 The rank-n simple EJA consisting of real symmetric n-by-n
1121 matrices, the usual symmetric Jordan product, and the trace inner
1122 product. It has dimension `(n^2 + n)/2` over the reals.
1126 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1130 sage: J = RealSymmetricEJA(2)
1131 sage: e0, e1, e2 = J.gens()
1139 In theory, our "field" can be any subfield of the reals::
1141 sage: RealSymmetricEJA(2, RDF)
1142 Euclidean Jordan algebra of dimension 3 over Real Double Field
1143 sage: RealSymmetricEJA(2, RR)
1144 Euclidean Jordan algebra of dimension 3 over Real Field with
1145 53 bits of precision
1149 The dimension of this algebra is `(n^2 + n) / 2`::
1151 sage: set_random_seed()
1152 sage: n_max = RealSymmetricEJA._max_test_case_size()
1153 sage: n = ZZ.random_element(1, n_max)
1154 sage: J = RealSymmetricEJA(n)
1155 sage: J.dimension() == (n^2 + n)/2
1158 The Jordan multiplication is what we think it is::
1160 sage: set_random_seed()
1161 sage: J = RealSymmetricEJA.random_instance()
1162 sage: x,y = J.random_elements(2)
1163 sage: actual = (x*y).natural_representation()
1164 sage: X = x.natural_representation()
1165 sage: Y = y.natural_representation()
1166 sage: expected = (X*Y + Y*X)/2
1167 sage: actual == expected
1169 sage: J(expected) == x*y
1172 We can change the generator prefix::
1174 sage: RealSymmetricEJA(3, prefix='q').gens()
1175 (q0, q1, q2, q3, q4, q5)
1177 Our natural basis is normalized with respect to the natural inner
1178 product unless we specify otherwise::
1180 sage: set_random_seed()
1181 sage: J = RealSymmetricEJA.random_instance()
1182 sage: all( b.norm() == 1 for b in J.gens() )
1185 Since our natural basis is normalized with respect to the natural
1186 inner product, and since we know that this algebra is an EJA, any
1187 left-multiplication operator's matrix will be symmetric because
1188 natural->EJA basis representation is an isometry and within the EJA
1189 the operator is self-adjoint by the Jordan axiom::
1191 sage: set_random_seed()
1192 sage: x = RealSymmetricEJA.random_instance().random_element()
1193 sage: x.operator().matrix().is_symmetric()
1196 We can construct the (trivial) algebra of rank zero::
1198 sage: RealSymmetricEJA(0)
1199 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1203 def _denormalized_basis(cls
, n
, field
):
1205 Return a basis for the space of real symmetric n-by-n matrices.
1209 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1213 sage: set_random_seed()
1214 sage: n = ZZ.random_element(1,5)
1215 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1216 sage: all( M.is_symmetric() for M in B)
1220 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1224 for j
in range(i
+1):
1225 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1229 Sij
= Eij
+ Eij
.transpose()
1235 def _max_test_case_size():
1236 return 4 # Dimension 10
1239 def __init__(self
, n
, field
=AA
, **kwargs
):
1240 basis
= self
._denormalized
_basis
(n
, field
)
1241 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, **kwargs
)
1242 self
.rank
.set_cache(n
)
1245 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1249 Embed the n-by-n complex matrix ``M`` into the space of real
1250 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1251 bi` to the block matrix ``[[a,b],[-b,a]]``.
1255 sage: from mjo.eja.eja_algebra import \
1256 ....: ComplexMatrixEuclideanJordanAlgebra
1260 sage: F = QuadraticField(-1, 'I')
1261 sage: x1 = F(4 - 2*i)
1262 sage: x2 = F(1 + 2*i)
1265 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1266 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1275 Embedding is a homomorphism (isomorphism, in fact)::
1277 sage: set_random_seed()
1278 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1279 sage: n = ZZ.random_element(n_max)
1280 sage: F = QuadraticField(-1, 'I')
1281 sage: X = random_matrix(F, n)
1282 sage: Y = random_matrix(F, n)
1283 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1284 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1285 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1292 raise ValueError("the matrix 'M' must be square")
1294 # We don't need any adjoined elements...
1295 field
= M
.base_ring().base_ring()
1299 a
= z
.list()[0] # real part, I guess
1300 b
= z
.list()[1] # imag part, I guess
1301 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1303 return matrix
.block(field
, n
, blocks
)
1307 def real_unembed(M
):
1309 The inverse of _embed_complex_matrix().
1313 sage: from mjo.eja.eja_algebra import \
1314 ....: ComplexMatrixEuclideanJordanAlgebra
1318 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1319 ....: [-2, 1, -4, 3],
1320 ....: [ 9, 10, 11, 12],
1321 ....: [-10, 9, -12, 11] ])
1322 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1324 [ 10*I + 9 12*I + 11]
1328 Unembedding is the inverse of embedding::
1330 sage: set_random_seed()
1331 sage: F = QuadraticField(-1, 'I')
1332 sage: M = random_matrix(F, 3)
1333 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1334 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1340 raise ValueError("the matrix 'M' must be square")
1341 if not n
.mod(2).is_zero():
1342 raise ValueError("the matrix 'M' must be a complex embedding")
1344 # If "M" was normalized, its base ring might have roots
1345 # adjoined and they can stick around after unembedding.
1346 field
= M
.base_ring()
1347 R
= PolynomialRing(field
, 'z')
1350 # Sage doesn't know how to embed AA into QQbar, i.e. how
1351 # to adjoin sqrt(-1) to AA.
1354 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1357 # Go top-left to bottom-right (reading order), converting every
1358 # 2-by-2 block we see to a single complex element.
1360 for k
in range(n
/2):
1361 for j
in range(n
/2):
1362 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1363 if submat
[0,0] != submat
[1,1]:
1364 raise ValueError('bad on-diagonal submatrix')
1365 if submat
[0,1] != -submat
[1,0]:
1366 raise ValueError('bad off-diagonal submatrix')
1367 z
= submat
[0,0] + submat
[0,1]*i
1370 return matrix(F
, n
/2, elements
)
1374 def natural_inner_product(cls
,X
,Y
):
1376 Compute a natural inner product in this algebra directly from
1381 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1385 This gives the same answer as the slow, default method implemented
1386 in :class:`MatrixEuclideanJordanAlgebra`::
1388 sage: set_random_seed()
1389 sage: J = ComplexHermitianEJA.random_instance()
1390 sage: x,y = J.random_elements(2)
1391 sage: Xe = x.natural_representation()
1392 sage: Ye = y.natural_representation()
1393 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1394 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1395 sage: expected = (X*Y).trace().real()
1396 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1397 sage: actual == expected
1401 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1404 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1406 The rank-n simple EJA consisting of complex Hermitian n-by-n
1407 matrices over the real numbers, the usual symmetric Jordan product,
1408 and the real-part-of-trace inner product. It has dimension `n^2` over
1413 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1417 In theory, our "field" can be any subfield of the reals::
1419 sage: ComplexHermitianEJA(2, RDF)
1420 Euclidean Jordan algebra of dimension 4 over Real Double Field
1421 sage: ComplexHermitianEJA(2, RR)
1422 Euclidean Jordan algebra of dimension 4 over Real Field with
1423 53 bits of precision
1427 The dimension of this algebra is `n^2`::
1429 sage: set_random_seed()
1430 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1431 sage: n = ZZ.random_element(1, n_max)
1432 sage: J = ComplexHermitianEJA(n)
1433 sage: J.dimension() == n^2
1436 The Jordan multiplication is what we think it is::
1438 sage: set_random_seed()
1439 sage: J = ComplexHermitianEJA.random_instance()
1440 sage: x,y = J.random_elements(2)
1441 sage: actual = (x*y).natural_representation()
1442 sage: X = x.natural_representation()
1443 sage: Y = y.natural_representation()
1444 sage: expected = (X*Y + Y*X)/2
1445 sage: actual == expected
1447 sage: J(expected) == x*y
1450 We can change the generator prefix::
1452 sage: ComplexHermitianEJA(2, prefix='z').gens()
1455 Our natural basis is normalized with respect to the natural inner
1456 product unless we specify otherwise::
1458 sage: set_random_seed()
1459 sage: J = ComplexHermitianEJA.random_instance()
1460 sage: all( b.norm() == 1 for b in J.gens() )
1463 Since our natural basis is normalized with respect to the natural
1464 inner product, and since we know that this algebra is an EJA, any
1465 left-multiplication operator's matrix will be symmetric because
1466 natural->EJA basis representation is an isometry and within the EJA
1467 the operator is self-adjoint by the Jordan axiom::
1469 sage: set_random_seed()
1470 sage: x = ComplexHermitianEJA.random_instance().random_element()
1471 sage: x.operator().matrix().is_symmetric()
1474 We can construct the (trivial) algebra of rank zero::
1476 sage: ComplexHermitianEJA(0)
1477 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1482 def _denormalized_basis(cls
, n
, field
):
1484 Returns a basis for the space of complex Hermitian n-by-n matrices.
1486 Why do we embed these? Basically, because all of numerical linear
1487 algebra assumes that you're working with vectors consisting of `n`
1488 entries from a field and scalars from the same field. There's no way
1489 to tell SageMath that (for example) the vectors contain complex
1490 numbers, while the scalar field is real.
1494 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1498 sage: set_random_seed()
1499 sage: n = ZZ.random_element(1,5)
1500 sage: field = QuadraticField(2, 'sqrt2')
1501 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1502 sage: all( M.is_symmetric() for M in B)
1506 R
= PolynomialRing(field
, 'z')
1508 F
= field
.extension(z
**2 + 1, 'I')
1511 # This is like the symmetric case, but we need to be careful:
1513 # * We want conjugate-symmetry, not just symmetry.
1514 # * The diagonal will (as a result) be real.
1518 for j
in range(i
+1):
1519 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1521 Sij
= cls
.real_embed(Eij
)
1524 # The second one has a minus because it's conjugated.
1525 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1527 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1530 # Since we embedded these, we can drop back to the "field" that we
1531 # started with instead of the complex extension "F".
1532 return ( s
.change_ring(field
) for s
in S
)
1535 def __init__(self
, n
, field
=AA
, **kwargs
):
1536 basis
= self
._denormalized
_basis
(n
,field
)
1537 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1538 self
.rank
.set_cache(n
)
1541 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1545 Embed the n-by-n quaternion matrix ``M`` into the space of real
1546 matrices of size 4n-by-4n by first sending each quaternion entry `z
1547 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1548 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1553 sage: from mjo.eja.eja_algebra import \
1554 ....: QuaternionMatrixEuclideanJordanAlgebra
1558 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1559 sage: i,j,k = Q.gens()
1560 sage: x = 1 + 2*i + 3*j + 4*k
1561 sage: M = matrix(Q, 1, [[x]])
1562 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1568 Embedding is a homomorphism (isomorphism, in fact)::
1570 sage: set_random_seed()
1571 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1572 sage: n = ZZ.random_element(n_max)
1573 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1574 sage: X = random_matrix(Q, n)
1575 sage: Y = random_matrix(Q, n)
1576 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1577 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1578 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1583 quaternions
= M
.base_ring()
1586 raise ValueError("the matrix 'M' must be square")
1588 F
= QuadraticField(-1, 'I')
1593 t
= z
.coefficient_tuple()
1598 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1599 [-c
+ d
*i
, a
- b
*i
]])
1600 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1601 blocks
.append(realM
)
1603 # We should have real entries by now, so use the realest field
1604 # we've got for the return value.
1605 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1610 def real_unembed(M
):
1612 The inverse of _embed_quaternion_matrix().
1616 sage: from mjo.eja.eja_algebra import \
1617 ....: QuaternionMatrixEuclideanJordanAlgebra
1621 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1622 ....: [-2, 1, -4, 3],
1623 ....: [-3, 4, 1, -2],
1624 ....: [-4, -3, 2, 1]])
1625 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1626 [1 + 2*i + 3*j + 4*k]
1630 Unembedding is the inverse of embedding::
1632 sage: set_random_seed()
1633 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1634 sage: M = random_matrix(Q, 3)
1635 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1636 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1642 raise ValueError("the matrix 'M' must be square")
1643 if not n
.mod(4).is_zero():
1644 raise ValueError("the matrix 'M' must be a quaternion embedding")
1646 # Use the base ring of the matrix to ensure that its entries can be
1647 # multiplied by elements of the quaternion algebra.
1648 field
= M
.base_ring()
1649 Q
= QuaternionAlgebra(field
,-1,-1)
1652 # Go top-left to bottom-right (reading order), converting every
1653 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1656 for l
in range(n
/4):
1657 for m
in range(n
/4):
1658 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1659 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1660 if submat
[0,0] != submat
[1,1].conjugate():
1661 raise ValueError('bad on-diagonal submatrix')
1662 if submat
[0,1] != -submat
[1,0].conjugate():
1663 raise ValueError('bad off-diagonal submatrix')
1664 z
= submat
[0,0].real()
1665 z
+= submat
[0,0].imag()*i
1666 z
+= submat
[0,1].real()*j
1667 z
+= submat
[0,1].imag()*k
1670 return matrix(Q
, n
/4, elements
)
1674 def natural_inner_product(cls
,X
,Y
):
1676 Compute a natural inner product in this algebra directly from
1681 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1685 This gives the same answer as the slow, default method implemented
1686 in :class:`MatrixEuclideanJordanAlgebra`::
1688 sage: set_random_seed()
1689 sage: J = QuaternionHermitianEJA.random_instance()
1690 sage: x,y = J.random_elements(2)
1691 sage: Xe = x.natural_representation()
1692 sage: Ye = y.natural_representation()
1693 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1694 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1695 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1696 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1697 sage: actual == expected
1701 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1704 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1706 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1707 matrices, the usual symmetric Jordan product, and the
1708 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1713 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1717 In theory, our "field" can be any subfield of the reals::
1719 sage: QuaternionHermitianEJA(2, RDF)
1720 Euclidean Jordan algebra of dimension 6 over Real Double Field
1721 sage: QuaternionHermitianEJA(2, RR)
1722 Euclidean Jordan algebra of dimension 6 over Real Field with
1723 53 bits of precision
1727 The dimension of this algebra is `2*n^2 - n`::
1729 sage: set_random_seed()
1730 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1731 sage: n = ZZ.random_element(1, n_max)
1732 sage: J = QuaternionHermitianEJA(n)
1733 sage: J.dimension() == 2*(n^2) - n
1736 The Jordan multiplication is what we think it is::
1738 sage: set_random_seed()
1739 sage: J = QuaternionHermitianEJA.random_instance()
1740 sage: x,y = J.random_elements(2)
1741 sage: actual = (x*y).natural_representation()
1742 sage: X = x.natural_representation()
1743 sage: Y = y.natural_representation()
1744 sage: expected = (X*Y + Y*X)/2
1745 sage: actual == expected
1747 sage: J(expected) == x*y
1750 We can change the generator prefix::
1752 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1753 (a0, a1, a2, a3, a4, a5)
1755 Our natural basis is normalized with respect to the natural inner
1756 product unless we specify otherwise::
1758 sage: set_random_seed()
1759 sage: J = QuaternionHermitianEJA.random_instance()
1760 sage: all( b.norm() == 1 for b in J.gens() )
1763 Since our natural basis is normalized with respect to the natural
1764 inner product, and since we know that this algebra is an EJA, any
1765 left-multiplication operator's matrix will be symmetric because
1766 natural->EJA basis representation is an isometry and within the EJA
1767 the operator is self-adjoint by the Jordan axiom::
1769 sage: set_random_seed()
1770 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1771 sage: x.operator().matrix().is_symmetric()
1774 We can construct the (trivial) algebra of rank zero::
1776 sage: QuaternionHermitianEJA(0)
1777 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1781 def _denormalized_basis(cls
, n
, field
):
1783 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1785 Why do we embed these? Basically, because all of numerical
1786 linear algebra assumes that you're working with vectors consisting
1787 of `n` entries from a field and scalars from the same field. There's
1788 no way to tell SageMath that (for example) the vectors contain
1789 complex numbers, while the scalar field is real.
1793 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1797 sage: set_random_seed()
1798 sage: n = ZZ.random_element(1,5)
1799 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1800 sage: all( M.is_symmetric() for M in B )
1804 Q
= QuaternionAlgebra(QQ
,-1,-1)
1807 # This is like the symmetric case, but we need to be careful:
1809 # * We want conjugate-symmetry, not just symmetry.
1810 # * The diagonal will (as a result) be real.
1814 for j
in range(i
+1):
1815 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1817 Sij
= cls
.real_embed(Eij
)
1820 # The second, third, and fourth ones have a minus
1821 # because they're conjugated.
1822 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1824 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1826 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1828 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1831 # Since we embedded these, we can drop back to the "field" that we
1832 # started with instead of the quaternion algebra "Q".
1833 return ( s
.change_ring(field
) for s
in S
)
1836 def __init__(self
, n
, field
=AA
, **kwargs
):
1837 basis
= self
._denormalized
_basis
(n
,field
)
1838 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1839 self
.rank
.set_cache(n
)
1842 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1844 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1845 with the half-trace inner product and jordan product ``x*y =
1846 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
1847 symmetric positive-definite "bilinear form" matrix. It has
1848 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
1849 when ``B`` is the identity matrix of order ``n-1``.
1853 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1854 ....: JordanSpinEJA)
1858 When no bilinear form is specified, the identity matrix is used,
1859 and the resulting algebra is the Jordan spin algebra::
1861 sage: J0 = BilinearFormEJA(3)
1862 sage: J1 = JordanSpinEJA(3)
1863 sage: J0.multiplication_table() == J0.multiplication_table()
1868 We can create a zero-dimensional algebra::
1870 sage: J = BilinearFormEJA(0)
1874 We can check the multiplication condition given in the Jordan, von
1875 Neumann, and Wigner paper (and also discussed on my "On the
1876 symmetry..." paper). Note that this relies heavily on the standard
1877 choice of basis, as does anything utilizing the bilinear form matrix::
1879 sage: set_random_seed()
1880 sage: n = ZZ.random_element(5)
1881 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
1882 sage: B = M.transpose()*M
1883 sage: J = BilinearFormEJA(n, B=B)
1884 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
1885 sage: V = J.vector_space()
1886 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
1887 ....: for ei in eis ]
1888 sage: actual = [ sis[i]*sis[j]
1889 ....: for i in range(n-1)
1890 ....: for j in range(n-1) ]
1891 sage: expected = [ J.one() if i == j else J.zero()
1892 ....: for i in range(n-1)
1893 ....: for j in range(n-1) ]
1894 sage: actual == expected
1897 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
1899 self
._B
= matrix
.identity(field
, max(0,n
-1))
1903 V
= VectorSpace(field
, n
)
1904 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1913 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
1914 zbar
= y0
*xbar
+ x0
*ybar
1915 z
= V([z0
] + zbar
.list())
1916 mult_table
[i
][j
] = z
1918 # The rank of this algebra is two, unless we're in a
1919 # one-dimensional ambient space (because the rank is bounded
1920 # by the ambient dimension).
1921 fdeja
= super(BilinearFormEJA
, self
)
1922 fdeja
.__init
__(field
, mult_table
, **kwargs
)
1923 self
.rank
.set_cache(min(n
,2))
1925 def inner_product(self
, x
, y
):
1927 Half of the trace inner product.
1929 This is defined so that the special case of the Jordan spin
1930 algebra gets the usual inner product.
1934 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1938 Ensure that this is one-half of the trace inner-product when
1939 the algebra isn't just the reals (when ``n`` isn't one). This
1940 is in Faraut and Koranyi, and also my "On the symmetry..."
1943 sage: set_random_seed()
1944 sage: n = ZZ.random_element(2,5)
1945 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
1946 sage: B = M.transpose()*M
1947 sage: J = BilinearFormEJA(n, B=B)
1948 sage: x = J.random_element()
1949 sage: y = J.random_element()
1950 sage: x.inner_product(y) == (x*y).trace()/2
1954 xvec
= x
.to_vector()
1956 yvec
= y
.to_vector()
1958 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
1961 class JordanSpinEJA(BilinearFormEJA
):
1963 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1964 with the usual inner product and jordan product ``x*y =
1965 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1970 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1974 This multiplication table can be verified by hand::
1976 sage: J = JordanSpinEJA(4)
1977 sage: e0,e1,e2,e3 = J.gens()
1993 We can change the generator prefix::
1995 sage: JordanSpinEJA(2, prefix='B').gens()
2000 Ensure that we have the usual inner product on `R^n`::
2002 sage: set_random_seed()
2003 sage: J = JordanSpinEJA.random_instance()
2004 sage: x,y = J.random_elements(2)
2005 sage: X = x.natural_representation()
2006 sage: Y = y.natural_representation()
2007 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2011 def __init__(self
, n
, field
=AA
, **kwargs
):
2012 # This is a special case of the BilinearFormEJA with the identity
2013 # matrix as its bilinear form.
2014 return super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2017 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2019 The trivial Euclidean Jordan algebra consisting of only a zero element.
2023 sage: from mjo.eja.eja_algebra import TrivialEJA
2027 sage: J = TrivialEJA()
2034 sage: 7*J.one()*12*J.one()
2036 sage: J.one().inner_product(J.one())
2038 sage: J.one().norm()
2040 sage: J.one().subalgebra_generated_by()
2041 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2046 def __init__(self
, field
=AA
, **kwargs
):
2048 fdeja
= super(TrivialEJA
, self
)
2049 # The rank is zero using my definition, namely the dimension of the
2050 # largest subalgebra generated by any element.
2051 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2052 self
.rank
.set_cache(0)