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 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
443 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
445 return self
._natural
_basis
[0].matrix_space()
449 def natural_inner_product(X
,Y
):
451 Compute the inner product of two naturally-represented elements.
453 For example in the real symmetric matrix EJA, this will compute
454 the trace inner-product of two n-by-n symmetric matrices. The
455 default should work for the real cartesian product EJA, the
456 Jordan spin EJA, and the real symmetric matrices. The others
457 will have to be overridden.
459 return (X
.conjugate_transpose()*Y
).trace()
465 Return the unit element of this algebra.
469 sage: from mjo.eja.eja_algebra import (HadamardEJA,
474 sage: J = HadamardEJA(5)
476 e0 + e1 + e2 + e3 + e4
480 The identity element acts like the identity::
482 sage: set_random_seed()
483 sage: J = random_eja()
484 sage: x = J.random_element()
485 sage: J.one()*x == x and x*J.one() == x
488 The matrix of the unit element's operator is the identity::
490 sage: set_random_seed()
491 sage: J = random_eja()
492 sage: actual = J.one().operator().matrix()
493 sage: expected = matrix.identity(J.base_ring(), J.dimension())
494 sage: actual == expected
498 # We can brute-force compute the matrices of the operators
499 # that correspond to the basis elements of this algebra.
500 # If some linear combination of those basis elements is the
501 # algebra identity, then the same linear combination of
502 # their matrices has to be the identity matrix.
504 # Of course, matrices aren't vectors in sage, so we have to
505 # appeal to the "long vectors" isometry.
506 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
508 # Now we use basis linear algebra to find the coefficients,
509 # of the matrices-as-vectors-linear-combination, which should
510 # work for the original algebra basis too.
511 A
= matrix
.column(self
.base_ring(), oper_vecs
)
513 # We used the isometry on the left-hand side already, but we
514 # still need to do it for the right-hand side. Recall that we
515 # wanted something that summed to the identity matrix.
516 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
518 # Now if there's an identity element in the algebra, this should work.
519 coeffs
= A
.solve_right(b
)
520 return self
.linear_combination(zip(self
.gens(), coeffs
))
523 def peirce_decomposition(self
, c
):
525 The Peirce decomposition of this algebra relative to the
528 In the future, this can be extended to a complete system of
529 orthogonal idempotents.
533 - ``c`` -- an idempotent of this algebra.
537 A triple (J0, J5, J1) containing two subalgebras and one subspace
540 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
541 corresponding to the eigenvalue zero.
543 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
544 corresponding to the eigenvalue one-half.
546 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
547 corresponding to the eigenvalue one.
549 These are the only possible eigenspaces for that operator, and this
550 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
551 orthogonal, and are subalgebras of this algebra with the appropriate
556 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
560 The canonical example comes from the symmetric matrices, which
561 decompose into diagonal and off-diagonal parts::
563 sage: J = RealSymmetricEJA(3)
564 sage: C = matrix(QQ, [ [1,0,0],
568 sage: J0,J5,J1 = J.peirce_decomposition(c)
570 Euclidean Jordan algebra of dimension 1...
572 Vector space of degree 6 and dimension 2...
574 Euclidean Jordan algebra of dimension 3...
578 Every algebra decomposes trivially with respect to its identity
581 sage: set_random_seed()
582 sage: J = random_eja()
583 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
584 sage: J0.dimension() == 0 and J5.dimension() == 0
586 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
589 The identity elements in the two subalgebras are the
590 projections onto their respective subspaces of the
591 superalgebra's identity element::
593 sage: set_random_seed()
594 sage: J = random_eja()
595 sage: x = J.random_element()
596 sage: if not J.is_trivial():
597 ....: while x.is_nilpotent():
598 ....: x = J.random_element()
599 sage: c = x.subalgebra_idempotent()
600 sage: J0,J5,J1 = J.peirce_decomposition(c)
601 sage: J1(c) == J1.one()
603 sage: J0(J.one() - c) == J0.one()
607 if not c
.is_idempotent():
608 raise ValueError("element is not idempotent: %s" % c
)
610 # Default these to what they should be if they turn out to be
611 # trivial, because eigenspaces_left() won't return eigenvalues
612 # corresponding to trivial spaces (e.g. it returns only the
613 # eigenspace corresponding to lambda=1 if you take the
614 # decomposition relative to the identity element).
615 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
616 J0
= trivial
# eigenvalue zero
617 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
618 J1
= trivial
# eigenvalue one
620 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
621 if eigval
== ~
(self
.base_ring()(2)):
624 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
625 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
631 raise ValueError("unexpected eigenvalue: %s" % eigval
)
636 def random_elements(self
, count
):
638 Return ``count`` random elements as a tuple.
642 sage: from mjo.eja.eja_algebra import JordanSpinEJA
646 sage: J = JordanSpinEJA(3)
647 sage: x,y,z = J.random_elements(3)
648 sage: all( [ x in J, y in J, z in J ])
650 sage: len( J.random_elements(10) ) == 10
654 return tuple( self
.random_element() for idx
in range(count
) )
657 def random_instance(cls
, field
=AA
, **kwargs
):
659 Return a random instance of this type of algebra.
661 Beware, this will crash for "most instances" because the
662 constructor below looks wrong.
664 if cls
is TrivialEJA
:
665 # The TrivialEJA class doesn't take an "n" argument because
669 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
670 return cls(n
, field
, **kwargs
)
673 def _charpoly_coefficients(self
):
675 The `r` polynomial coefficients of the "characteristic polynomial
679 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
680 R
= PolynomialRing(self
.base_ring(), var_names
)
682 F
= R
.fraction_field()
685 # From a result in my book, these are the entries of the
686 # basis representation of L_x.
687 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
690 L_x
= matrix(F
, n
, n
, L_x_i_j
)
693 if self
.rank
.is_in_cache():
695 # There's no need to pad the system with redundant
696 # columns if we *know* they'll be redundant.
699 # Compute an extra power in case the rank is equal to
700 # the dimension (otherwise, we would stop at x^(r-1)).
701 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
702 for k
in range(n
+1) ]
703 A
= matrix
.column(F
, x_powers
[:n
])
704 AE
= A
.extended_echelon_form()
711 # The theory says that only the first "r" coefficients are
712 # nonzero, and they actually live in the original polynomial
713 # ring and not the fraction field. We negate them because
714 # in the actual characteristic polynomial, they get moved
715 # to the other side where x^r lives.
716 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
721 Return the rank of this EJA.
723 This is a cached method because we know the rank a priori for
724 all of the algebras we can construct. Thus we can avoid the
725 expensive ``_charpoly_coefficients()`` call unless we truly
726 need to compute the whole characteristic polynomial.
730 sage: from mjo.eja.eja_algebra import (HadamardEJA,
732 ....: RealSymmetricEJA,
733 ....: ComplexHermitianEJA,
734 ....: QuaternionHermitianEJA,
739 The rank of the Jordan spin algebra is always two::
741 sage: JordanSpinEJA(2).rank()
743 sage: JordanSpinEJA(3).rank()
745 sage: JordanSpinEJA(4).rank()
748 The rank of the `n`-by-`n` Hermitian real, complex, or
749 quaternion matrices is `n`::
751 sage: RealSymmetricEJA(4).rank()
753 sage: ComplexHermitianEJA(3).rank()
755 sage: QuaternionHermitianEJA(2).rank()
760 Ensure that every EJA that we know how to construct has a
761 positive integer rank, unless the algebra is trivial in
762 which case its rank will be zero::
764 sage: set_random_seed()
765 sage: J = random_eja()
769 sage: r > 0 or (r == 0 and J.is_trivial())
772 Ensure that computing the rank actually works, since the ranks
773 of all simple algebras are known and will be cached by default::
775 sage: J = HadamardEJA(4)
776 sage: J.rank.clear_cache()
782 sage: J = JordanSpinEJA(4)
783 sage: J.rank.clear_cache()
789 sage: J = RealSymmetricEJA(3)
790 sage: J.rank.clear_cache()
796 sage: J = ComplexHermitianEJA(2)
797 sage: J.rank.clear_cache()
803 sage: J = QuaternionHermitianEJA(2)
804 sage: J.rank.clear_cache()
808 return len(self
._charpoly
_coefficients
())
811 def vector_space(self
):
813 Return the vector space that underlies this algebra.
817 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
821 sage: J = RealSymmetricEJA(2)
822 sage: J.vector_space()
823 Vector space of dimension 3 over...
826 return self
.zero().to_vector().parent().ambient_vector_space()
829 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
832 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
):
834 Return the Euclidean Jordan Algebra corresponding to the set
835 `R^n` under the Hadamard product.
837 Note: this is nothing more than the Cartesian product of ``n``
838 copies of the spin algebra. Once Cartesian product algebras
839 are implemented, this can go.
843 sage: from mjo.eja.eja_algebra import HadamardEJA
847 This multiplication table can be verified by hand::
849 sage: J = HadamardEJA(3)
850 sage: e0,e1,e2 = J.gens()
866 We can change the generator prefix::
868 sage: HadamardEJA(3, prefix='r').gens()
872 def __init__(self
, n
, field
=AA
, **kwargs
):
873 V
= VectorSpace(field
, n
)
874 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
877 fdeja
= super(HadamardEJA
, self
)
878 fdeja
.__init
__(field
, mult_table
, **kwargs
)
879 self
.rank
.set_cache(n
)
881 def inner_product(self
, x
, y
):
883 Faster to reimplement than to use natural representations.
887 sage: from mjo.eja.eja_algebra import HadamardEJA
891 Ensure that this is the usual inner product for the algebras
894 sage: set_random_seed()
895 sage: J = HadamardEJA.random_instance()
896 sage: x,y = J.random_elements(2)
897 sage: X = x.natural_representation()
898 sage: Y = y.natural_representation()
899 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
903 return x
.to_vector().inner_product(y
.to_vector())
906 def random_eja(field
=AA
, nontrivial
=False):
908 Return a "random" finite-dimensional Euclidean Jordan Algebra.
912 sage: from mjo.eja.eja_algebra import random_eja
917 Euclidean Jordan algebra of dimension...
920 eja_classes
= [HadamardEJA
,
924 QuaternionHermitianEJA
]
926 eja_classes
.append(TrivialEJA
)
927 classname
= choice(eja_classes
)
928 return classname
.random_instance(field
=field
)
935 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
937 def _max_test_case_size():
938 # Play it safe, since this will be squared and the underlying
939 # field can have dimension 4 (quaternions) too.
942 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
944 Compared to the superclass constructor, we take a basis instead of
945 a multiplication table because the latter can be computed in terms
946 of the former when the product is known (like it is here).
948 # Used in this class's fast _charpoly_coefficients() override.
949 self
._basis
_normalizers
= None
951 # We're going to loop through this a few times, so now's a good
952 # time to ensure that it isn't a generator expression.
955 if len(basis
) > 1 and normalize_basis
:
956 # We'll need sqrt(2) to normalize the basis, and this
957 # winds up in the multiplication table, so the whole
958 # algebra needs to be over the field extension.
959 R
= PolynomialRing(field
, 'z')
962 if p
.is_irreducible():
963 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
964 basis
= tuple( s
.change_ring(field
) for s
in basis
)
965 self
._basis
_normalizers
= tuple(
966 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
967 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
969 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
971 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
972 fdeja
.__init
__(field
, Qs
, natural_basis
=basis
, **kwargs
)
977 def _charpoly_coefficients(self
):
979 Override the parent method with something that tries to compute
980 over a faster (non-extension) field.
982 if self
._basis
_normalizers
is None:
983 # We didn't normalize, so assume that the basis we started
984 # with had entries in a nice field.
985 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
987 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
988 self
._basis
_normalizers
) )
990 # Do this over the rationals and convert back at the end.
991 # Only works because we know the entries of the basis are
993 J
= MatrixEuclideanJordanAlgebra(QQ
,
995 normalize_basis
=False)
996 a
= J
._charpoly
_coefficients
()
998 # Unfortunately, changing the basis does change the
999 # coefficients of the characteristic polynomial, but since
1000 # these are really the coefficients of the "characteristic
1001 # polynomial of" function, everything is still nice and
1002 # unevaluated. It's therefore "obvious" how scaling the
1003 # basis affects the coordinate variables X1, X2, et
1004 # cetera. Scaling the first basis vector up by "n" adds a
1005 # factor of 1/n into every "X1" term, for example. So here
1006 # we simply undo the basis_normalizer scaling that we
1007 # performed earlier.
1009 # The a[0] access here is safe because trivial algebras
1010 # won't have any basis normalizers and therefore won't
1011 # make it to this "else" branch.
1012 XS
= a
[0].parent().gens()
1013 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1014 for i
in range(len(XS
)) }
1015 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1019 def multiplication_table_from_matrix_basis(basis
):
1021 At least three of the five simple Euclidean Jordan algebras have the
1022 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1023 multiplication on the right is matrix multiplication. Given a basis
1024 for the underlying matrix space, this function returns a
1025 multiplication table (obtained by looping through the basis
1026 elements) for an algebra of those matrices.
1028 # In S^2, for example, we nominally have four coordinates even
1029 # though the space is of dimension three only. The vector space V
1030 # is supposed to hold the entire long vector, and the subspace W
1031 # of V will be spanned by the vectors that arise from symmetric
1032 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1036 field
= basis
[0].base_ring()
1037 dimension
= basis
[0].nrows()
1039 V
= VectorSpace(field
, dimension
**2)
1040 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1042 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1045 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1046 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1054 Embed the matrix ``M`` into a space of real matrices.
1056 The matrix ``M`` can have entries in any field at the moment:
1057 the real numbers, complex numbers, or quaternions. And although
1058 they are not a field, we can probably support octonions at some
1059 point, too. This function returns a real matrix that "acts like"
1060 the original with respect to matrix multiplication; i.e.
1062 real_embed(M*N) = real_embed(M)*real_embed(N)
1065 raise NotImplementedError
1069 def real_unembed(M
):
1071 The inverse of :meth:`real_embed`.
1073 raise NotImplementedError
1077 def natural_inner_product(cls
,X
,Y
):
1078 Xu
= cls
.real_unembed(X
)
1079 Yu
= cls
.real_unembed(Y
)
1080 tr
= (Xu
*Yu
).trace()
1083 # It's real already.
1086 # Otherwise, try the thing that works for complex numbers; and
1087 # if that doesn't work, the thing that works for quaternions.
1089 return tr
.vector()[0] # real part, imag part is index 1
1090 except AttributeError:
1091 # A quaternions doesn't have a vector() method, but does
1092 # have coefficient_tuple() method that returns the
1093 # coefficients of 1, i, j, and k -- in that order.
1094 return tr
.coefficient_tuple()[0]
1097 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1101 The identity function, for embedding real matrices into real
1107 def real_unembed(M
):
1109 The identity function, for unembedding real matrices from real
1115 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1117 The rank-n simple EJA consisting of real symmetric n-by-n
1118 matrices, the usual symmetric Jordan product, and the trace inner
1119 product. It has dimension `(n^2 + n)/2` over the reals.
1123 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1127 sage: J = RealSymmetricEJA(2)
1128 sage: e0, e1, e2 = J.gens()
1136 In theory, our "field" can be any subfield of the reals::
1138 sage: RealSymmetricEJA(2, RDF)
1139 Euclidean Jordan algebra of dimension 3 over Real Double Field
1140 sage: RealSymmetricEJA(2, RR)
1141 Euclidean Jordan algebra of dimension 3 over Real Field with
1142 53 bits of precision
1146 The dimension of this algebra is `(n^2 + n) / 2`::
1148 sage: set_random_seed()
1149 sage: n_max = RealSymmetricEJA._max_test_case_size()
1150 sage: n = ZZ.random_element(1, n_max)
1151 sage: J = RealSymmetricEJA(n)
1152 sage: J.dimension() == (n^2 + n)/2
1155 The Jordan multiplication is what we think it is::
1157 sage: set_random_seed()
1158 sage: J = RealSymmetricEJA.random_instance()
1159 sage: x,y = J.random_elements(2)
1160 sage: actual = (x*y).natural_representation()
1161 sage: X = x.natural_representation()
1162 sage: Y = y.natural_representation()
1163 sage: expected = (X*Y + Y*X)/2
1164 sage: actual == expected
1166 sage: J(expected) == x*y
1169 We can change the generator prefix::
1171 sage: RealSymmetricEJA(3, prefix='q').gens()
1172 (q0, q1, q2, q3, q4, q5)
1174 Our natural basis is normalized with respect to the natural inner
1175 product unless we specify otherwise::
1177 sage: set_random_seed()
1178 sage: J = RealSymmetricEJA.random_instance()
1179 sage: all( b.norm() == 1 for b in J.gens() )
1182 Since our natural basis is normalized with respect to the natural
1183 inner product, and since we know that this algebra is an EJA, any
1184 left-multiplication operator's matrix will be symmetric because
1185 natural->EJA basis representation is an isometry and within the EJA
1186 the operator is self-adjoint by the Jordan axiom::
1188 sage: set_random_seed()
1189 sage: x = RealSymmetricEJA.random_instance().random_element()
1190 sage: x.operator().matrix().is_symmetric()
1193 We can construct the (trivial) algebra of rank zero::
1195 sage: RealSymmetricEJA(0)
1196 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1200 def _denormalized_basis(cls
, n
, field
):
1202 Return a basis for the space of real symmetric n-by-n matrices.
1206 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1210 sage: set_random_seed()
1211 sage: n = ZZ.random_element(1,5)
1212 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1213 sage: all( M.is_symmetric() for M in B)
1217 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1221 for j
in range(i
+1):
1222 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1226 Sij
= Eij
+ Eij
.transpose()
1232 def _max_test_case_size():
1233 return 4 # Dimension 10
1236 def __init__(self
, n
, field
=AA
, **kwargs
):
1237 basis
= self
._denormalized
_basis
(n
, field
)
1238 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, **kwargs
)
1239 self
.rank
.set_cache(n
)
1242 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1246 Embed the n-by-n complex matrix ``M`` into the space of real
1247 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1248 bi` to the block matrix ``[[a,b],[-b,a]]``.
1252 sage: from mjo.eja.eja_algebra import \
1253 ....: ComplexMatrixEuclideanJordanAlgebra
1257 sage: F = QuadraticField(-1, 'I')
1258 sage: x1 = F(4 - 2*i)
1259 sage: x2 = F(1 + 2*i)
1262 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1263 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1272 Embedding is a homomorphism (isomorphism, in fact)::
1274 sage: set_random_seed()
1275 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1276 sage: n = ZZ.random_element(n_max)
1277 sage: F = QuadraticField(-1, 'I')
1278 sage: X = random_matrix(F, n)
1279 sage: Y = random_matrix(F, n)
1280 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1281 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1282 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1289 raise ValueError("the matrix 'M' must be square")
1291 # We don't need any adjoined elements...
1292 field
= M
.base_ring().base_ring()
1296 a
= z
.list()[0] # real part, I guess
1297 b
= z
.list()[1] # imag part, I guess
1298 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1300 return matrix
.block(field
, n
, blocks
)
1304 def real_unembed(M
):
1306 The inverse of _embed_complex_matrix().
1310 sage: from mjo.eja.eja_algebra import \
1311 ....: ComplexMatrixEuclideanJordanAlgebra
1315 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1316 ....: [-2, 1, -4, 3],
1317 ....: [ 9, 10, 11, 12],
1318 ....: [-10, 9, -12, 11] ])
1319 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1321 [ 10*I + 9 12*I + 11]
1325 Unembedding is the inverse of embedding::
1327 sage: set_random_seed()
1328 sage: F = QuadraticField(-1, 'I')
1329 sage: M = random_matrix(F, 3)
1330 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1331 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1337 raise ValueError("the matrix 'M' must be square")
1338 if not n
.mod(2).is_zero():
1339 raise ValueError("the matrix 'M' must be a complex embedding")
1341 # If "M" was normalized, its base ring might have roots
1342 # adjoined and they can stick around after unembedding.
1343 field
= M
.base_ring()
1344 R
= PolynomialRing(field
, 'z')
1347 # Sage doesn't know how to embed AA into QQbar, i.e. how
1348 # to adjoin sqrt(-1) to AA.
1351 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1354 # Go top-left to bottom-right (reading order), converting every
1355 # 2-by-2 block we see to a single complex element.
1357 for k
in range(n
/2):
1358 for j
in range(n
/2):
1359 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1360 if submat
[0,0] != submat
[1,1]:
1361 raise ValueError('bad on-diagonal submatrix')
1362 if submat
[0,1] != -submat
[1,0]:
1363 raise ValueError('bad off-diagonal submatrix')
1364 z
= submat
[0,0] + submat
[0,1]*i
1367 return matrix(F
, n
/2, elements
)
1371 def natural_inner_product(cls
,X
,Y
):
1373 Compute a natural inner product in this algebra directly from
1378 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1382 This gives the same answer as the slow, default method implemented
1383 in :class:`MatrixEuclideanJordanAlgebra`::
1385 sage: set_random_seed()
1386 sage: J = ComplexHermitianEJA.random_instance()
1387 sage: x,y = J.random_elements(2)
1388 sage: Xe = x.natural_representation()
1389 sage: Ye = y.natural_representation()
1390 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1391 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1392 sage: expected = (X*Y).trace().real()
1393 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1394 sage: actual == expected
1398 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1401 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1403 The rank-n simple EJA consisting of complex Hermitian n-by-n
1404 matrices over the real numbers, the usual symmetric Jordan product,
1405 and the real-part-of-trace inner product. It has dimension `n^2` over
1410 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1414 In theory, our "field" can be any subfield of the reals::
1416 sage: ComplexHermitianEJA(2, RDF)
1417 Euclidean Jordan algebra of dimension 4 over Real Double Field
1418 sage: ComplexHermitianEJA(2, RR)
1419 Euclidean Jordan algebra of dimension 4 over Real Field with
1420 53 bits of precision
1424 The dimension of this algebra is `n^2`::
1426 sage: set_random_seed()
1427 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1428 sage: n = ZZ.random_element(1, n_max)
1429 sage: J = ComplexHermitianEJA(n)
1430 sage: J.dimension() == n^2
1433 The Jordan multiplication is what we think it is::
1435 sage: set_random_seed()
1436 sage: J = ComplexHermitianEJA.random_instance()
1437 sage: x,y = J.random_elements(2)
1438 sage: actual = (x*y).natural_representation()
1439 sage: X = x.natural_representation()
1440 sage: Y = y.natural_representation()
1441 sage: expected = (X*Y + Y*X)/2
1442 sage: actual == expected
1444 sage: J(expected) == x*y
1447 We can change the generator prefix::
1449 sage: ComplexHermitianEJA(2, prefix='z').gens()
1452 Our natural basis is normalized with respect to the natural inner
1453 product unless we specify otherwise::
1455 sage: set_random_seed()
1456 sage: J = ComplexHermitianEJA.random_instance()
1457 sage: all( b.norm() == 1 for b in J.gens() )
1460 Since our natural basis is normalized with respect to the natural
1461 inner product, and since we know that this algebra is an EJA, any
1462 left-multiplication operator's matrix will be symmetric because
1463 natural->EJA basis representation is an isometry and within the EJA
1464 the operator is self-adjoint by the Jordan axiom::
1466 sage: set_random_seed()
1467 sage: x = ComplexHermitianEJA.random_instance().random_element()
1468 sage: x.operator().matrix().is_symmetric()
1471 We can construct the (trivial) algebra of rank zero::
1473 sage: ComplexHermitianEJA(0)
1474 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1479 def _denormalized_basis(cls
, n
, field
):
1481 Returns a basis for the space of complex Hermitian n-by-n matrices.
1483 Why do we embed these? Basically, because all of numerical linear
1484 algebra assumes that you're working with vectors consisting of `n`
1485 entries from a field and scalars from the same field. There's no way
1486 to tell SageMath that (for example) the vectors contain complex
1487 numbers, while the scalar field is real.
1491 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1495 sage: set_random_seed()
1496 sage: n = ZZ.random_element(1,5)
1497 sage: field = QuadraticField(2, 'sqrt2')
1498 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1499 sage: all( M.is_symmetric() for M in B)
1503 R
= PolynomialRing(field
, 'z')
1505 F
= field
.extension(z
**2 + 1, 'I')
1508 # This is like the symmetric case, but we need to be careful:
1510 # * We want conjugate-symmetry, not just symmetry.
1511 # * The diagonal will (as a result) be real.
1515 for j
in range(i
+1):
1516 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1518 Sij
= cls
.real_embed(Eij
)
1521 # The second one has a minus because it's conjugated.
1522 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1524 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1527 # Since we embedded these, we can drop back to the "field" that we
1528 # started with instead of the complex extension "F".
1529 return ( s
.change_ring(field
) for s
in S
)
1532 def __init__(self
, n
, field
=AA
, **kwargs
):
1533 basis
= self
._denormalized
_basis
(n
,field
)
1534 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1535 self
.rank
.set_cache(n
)
1538 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1542 Embed the n-by-n quaternion matrix ``M`` into the space of real
1543 matrices of size 4n-by-4n by first sending each quaternion entry `z
1544 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1545 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1550 sage: from mjo.eja.eja_algebra import \
1551 ....: QuaternionMatrixEuclideanJordanAlgebra
1555 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1556 sage: i,j,k = Q.gens()
1557 sage: x = 1 + 2*i + 3*j + 4*k
1558 sage: M = matrix(Q, 1, [[x]])
1559 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1565 Embedding is a homomorphism (isomorphism, in fact)::
1567 sage: set_random_seed()
1568 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1569 sage: n = ZZ.random_element(n_max)
1570 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1571 sage: X = random_matrix(Q, n)
1572 sage: Y = random_matrix(Q, n)
1573 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1574 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1575 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1580 quaternions
= M
.base_ring()
1583 raise ValueError("the matrix 'M' must be square")
1585 F
= QuadraticField(-1, 'I')
1590 t
= z
.coefficient_tuple()
1595 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1596 [-c
+ d
*i
, a
- b
*i
]])
1597 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1598 blocks
.append(realM
)
1600 # We should have real entries by now, so use the realest field
1601 # we've got for the return value.
1602 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1607 def real_unembed(M
):
1609 The inverse of _embed_quaternion_matrix().
1613 sage: from mjo.eja.eja_algebra import \
1614 ....: QuaternionMatrixEuclideanJordanAlgebra
1618 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1619 ....: [-2, 1, -4, 3],
1620 ....: [-3, 4, 1, -2],
1621 ....: [-4, -3, 2, 1]])
1622 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1623 [1 + 2*i + 3*j + 4*k]
1627 Unembedding is the inverse of embedding::
1629 sage: set_random_seed()
1630 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1631 sage: M = random_matrix(Q, 3)
1632 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1633 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1639 raise ValueError("the matrix 'M' must be square")
1640 if not n
.mod(4).is_zero():
1641 raise ValueError("the matrix 'M' must be a quaternion embedding")
1643 # Use the base ring of the matrix to ensure that its entries can be
1644 # multiplied by elements of the quaternion algebra.
1645 field
= M
.base_ring()
1646 Q
= QuaternionAlgebra(field
,-1,-1)
1649 # Go top-left to bottom-right (reading order), converting every
1650 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1653 for l
in range(n
/4):
1654 for m
in range(n
/4):
1655 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1656 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1657 if submat
[0,0] != submat
[1,1].conjugate():
1658 raise ValueError('bad on-diagonal submatrix')
1659 if submat
[0,1] != -submat
[1,0].conjugate():
1660 raise ValueError('bad off-diagonal submatrix')
1661 z
= submat
[0,0].real()
1662 z
+= submat
[0,0].imag()*i
1663 z
+= submat
[0,1].real()*j
1664 z
+= submat
[0,1].imag()*k
1667 return matrix(Q
, n
/4, elements
)
1671 def natural_inner_product(cls
,X
,Y
):
1673 Compute a natural inner product in this algebra directly from
1678 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1682 This gives the same answer as the slow, default method implemented
1683 in :class:`MatrixEuclideanJordanAlgebra`::
1685 sage: set_random_seed()
1686 sage: J = QuaternionHermitianEJA.random_instance()
1687 sage: x,y = J.random_elements(2)
1688 sage: Xe = x.natural_representation()
1689 sage: Ye = y.natural_representation()
1690 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1691 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1692 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1693 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1694 sage: actual == expected
1698 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1701 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1703 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1704 matrices, the usual symmetric Jordan product, and the
1705 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1710 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1714 In theory, our "field" can be any subfield of the reals::
1716 sage: QuaternionHermitianEJA(2, RDF)
1717 Euclidean Jordan algebra of dimension 6 over Real Double Field
1718 sage: QuaternionHermitianEJA(2, RR)
1719 Euclidean Jordan algebra of dimension 6 over Real Field with
1720 53 bits of precision
1724 The dimension of this algebra is `2*n^2 - n`::
1726 sage: set_random_seed()
1727 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1728 sage: n = ZZ.random_element(1, n_max)
1729 sage: J = QuaternionHermitianEJA(n)
1730 sage: J.dimension() == 2*(n^2) - n
1733 The Jordan multiplication is what we think it is::
1735 sage: set_random_seed()
1736 sage: J = QuaternionHermitianEJA.random_instance()
1737 sage: x,y = J.random_elements(2)
1738 sage: actual = (x*y).natural_representation()
1739 sage: X = x.natural_representation()
1740 sage: Y = y.natural_representation()
1741 sage: expected = (X*Y + Y*X)/2
1742 sage: actual == expected
1744 sage: J(expected) == x*y
1747 We can change the generator prefix::
1749 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1750 (a0, a1, a2, a3, a4, a5)
1752 Our natural basis is normalized with respect to the natural inner
1753 product unless we specify otherwise::
1755 sage: set_random_seed()
1756 sage: J = QuaternionHermitianEJA.random_instance()
1757 sage: all( b.norm() == 1 for b in J.gens() )
1760 Since our natural basis is normalized with respect to the natural
1761 inner product, and since we know that this algebra is an EJA, any
1762 left-multiplication operator's matrix will be symmetric because
1763 natural->EJA basis representation is an isometry and within the EJA
1764 the operator is self-adjoint by the Jordan axiom::
1766 sage: set_random_seed()
1767 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1768 sage: x.operator().matrix().is_symmetric()
1771 We can construct the (trivial) algebra of rank zero::
1773 sage: QuaternionHermitianEJA(0)
1774 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1778 def _denormalized_basis(cls
, n
, field
):
1780 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1782 Why do we embed these? Basically, because all of numerical
1783 linear algebra assumes that you're working with vectors consisting
1784 of `n` entries from a field and scalars from the same field. There's
1785 no way to tell SageMath that (for example) the vectors contain
1786 complex numbers, while the scalar field is real.
1790 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1794 sage: set_random_seed()
1795 sage: n = ZZ.random_element(1,5)
1796 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1797 sage: all( M.is_symmetric() for M in B )
1801 Q
= QuaternionAlgebra(QQ
,-1,-1)
1804 # This is like the symmetric case, but we need to be careful:
1806 # * We want conjugate-symmetry, not just symmetry.
1807 # * The diagonal will (as a result) be real.
1811 for j
in range(i
+1):
1812 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1814 Sij
= cls
.real_embed(Eij
)
1817 # The second, third, and fourth ones have a minus
1818 # because they're conjugated.
1819 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1821 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1823 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1825 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1828 # Since we embedded these, we can drop back to the "field" that we
1829 # started with instead of the quaternion algebra "Q".
1830 return ( s
.change_ring(field
) for s
in S
)
1833 def __init__(self
, n
, field
=AA
, **kwargs
):
1834 basis
= self
._denormalized
_basis
(n
,field
)
1835 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1836 self
.rank
.set_cache(n
)
1839 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1841 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1842 with the half-trace inner product and jordan product ``x*y =
1843 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
1844 symmetric positive-definite "bilinear form" matrix. It has
1845 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
1846 when ``B`` is the identity matrix of order ``n-1``.
1850 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1851 ....: JordanSpinEJA)
1855 When no bilinear form is specified, the identity matrix is used,
1856 and the resulting algebra is the Jordan spin algebra::
1858 sage: J0 = BilinearFormEJA(3)
1859 sage: J1 = JordanSpinEJA(3)
1860 sage: J0.multiplication_table() == J0.multiplication_table()
1865 We can create a zero-dimensional algebra::
1867 sage: J = BilinearFormEJA(0)
1871 We can check the multiplication condition given in the Jordan, von
1872 Neumann, and Wigner paper (and also discussed on my "On the
1873 symmetry..." paper). Note that this relies heavily on the standard
1874 choice of basis, as does anything utilizing the bilinear form matrix::
1876 sage: set_random_seed()
1877 sage: n = ZZ.random_element(5)
1878 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
1879 sage: B = M.transpose()*M
1880 sage: J = BilinearFormEJA(n, B=B)
1881 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
1882 sage: V = J.vector_space()
1883 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
1884 ....: for ei in eis ]
1885 sage: actual = [ sis[i]*sis[j]
1886 ....: for i in range(n-1)
1887 ....: for j in range(n-1) ]
1888 sage: expected = [ J.one() if i == j else J.zero()
1889 ....: for i in range(n-1)
1890 ....: for j in range(n-1) ]
1891 sage: actual == expected
1894 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
1896 self
._B
= matrix
.identity(field
, max(0,n
-1))
1900 V
= VectorSpace(field
, n
)
1901 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1910 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
1911 zbar
= y0
*xbar
+ x0
*ybar
1912 z
= V([z0
] + zbar
.list())
1913 mult_table
[i
][j
] = z
1915 # The rank of this algebra is two, unless we're in a
1916 # one-dimensional ambient space (because the rank is bounded
1917 # by the ambient dimension).
1918 fdeja
= super(BilinearFormEJA
, self
)
1919 fdeja
.__init
__(field
, mult_table
, **kwargs
)
1920 self
.rank
.set_cache(min(n
,2))
1922 def inner_product(self
, x
, y
):
1924 Half of the trace inner product.
1926 This is defined so that the special case of the Jordan spin
1927 algebra gets the usual inner product.
1931 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1935 Ensure that this is one-half of the trace inner-product when
1936 the algebra isn't just the reals (when ``n`` isn't one). This
1937 is in Faraut and Koranyi, and also my "On the symmetry..."
1940 sage: set_random_seed()
1941 sage: n = ZZ.random_element(2,5)
1942 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
1943 sage: B = M.transpose()*M
1944 sage: J = BilinearFormEJA(n, B=B)
1945 sage: x = J.random_element()
1946 sage: y = J.random_element()
1947 sage: x.inner_product(y) == (x*y).trace()/2
1951 xvec
= x
.to_vector()
1953 yvec
= y
.to_vector()
1955 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
1958 class JordanSpinEJA(BilinearFormEJA
):
1960 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1961 with the usual inner product and jordan product ``x*y =
1962 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1967 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1971 This multiplication table can be verified by hand::
1973 sage: J = JordanSpinEJA(4)
1974 sage: e0,e1,e2,e3 = J.gens()
1990 We can change the generator prefix::
1992 sage: JordanSpinEJA(2, prefix='B').gens()
1997 Ensure that we have the usual inner product on `R^n`::
1999 sage: set_random_seed()
2000 sage: J = JordanSpinEJA.random_instance()
2001 sage: x,y = J.random_elements(2)
2002 sage: X = x.natural_representation()
2003 sage: Y = y.natural_representation()
2004 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2008 def __init__(self
, n
, field
=AA
, **kwargs
):
2009 # This is a special case of the BilinearFormEJA with the identity
2010 # matrix as its bilinear form.
2011 return super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2014 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2016 The trivial Euclidean Jordan algebra consisting of only a zero element.
2020 sage: from mjo.eja.eja_algebra import TrivialEJA
2024 sage: J = TrivialEJA()
2031 sage: 7*J.one()*12*J.one()
2033 sage: J.one().inner_product(J.one())
2035 sage: J.one().norm()
2037 sage: J.one().subalgebra_generated_by()
2038 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2043 def __init__(self
, field
=AA
, **kwargs
):
2045 fdeja
= super(TrivialEJA
, self
)
2046 # The rank is zero using my definition, namely the dimension of the
2047 # largest subalgebra generated by any element.
2048 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2049 self
.rank
.set_cache(0)