2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 from itertools
import repeat
10 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
11 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
12 from sage
.combinat
.free_module
import CombinatorialFreeModule
13 from sage
.matrix
.constructor
import matrix
14 from sage
.matrix
.matrix_space
import MatrixSpace
15 from sage
.misc
.cachefunc
import cached_method
16 from sage
.misc
.lazy_import
import lazy_import
17 from sage
.misc
.prandom
import choice
18 from sage
.misc
.table
import table
19 from sage
.modules
.free_module
import FreeModule
, VectorSpace
20 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
23 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
24 lazy_import('mjo.eja.eja_subalgebra',
25 'FiniteDimensionalEuclideanJordanSubalgebra')
26 from mjo
.eja
.eja_utils
import _mat2vec
28 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
30 def _coerce_map_from_base_ring(self
):
32 Disable the map from the base ring into the algebra.
34 Performing a nonsense conversion like this automatically
35 is counterpedagogical. The fallback is to try the usual
36 element constructor, which should also fail.
40 sage: from mjo.eja.eja_algebra import random_eja
44 sage: set_random_seed()
45 sage: J = random_eja()
47 Traceback (most recent call last):
49 ValueError: not a naturally-represented algebra element
64 sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
68 By definition, Jordan multiplication commutes::
70 sage: set_random_seed()
71 sage: J = random_eja()
72 sage: x,y = J.random_elements(2)
78 The ``field`` we're given must be real::
80 sage: JordanSpinEJA(2,QQbar)
81 Traceback (most recent call last):
83 ValueError: field is not real
87 if not field
.is_subring(RR
):
88 # Note: this does return true for the real algebraic
89 # field, and any quadratic field where we've specified
91 raise ValueError('field is not real')
93 self
._natural
_basis
= natural_basis
96 category
= MagmaticAlgebras(field
).FiniteDimensional()
97 category
= category
.WithBasis().Unital()
99 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
101 range(len(mult_table
)),
104 self
.print_options(bracket
='')
106 # The multiplication table we're given is necessarily in terms
107 # of vectors, because we don't have an algebra yet for
108 # anything to be an element of. However, it's faster in the
109 # long run to have the multiplication table be in terms of
110 # algebra elements. We do this after calling the superclass
111 # constructor so that from_vector() knows what to do.
112 self
._multiplication
_table
= [
113 list(map(lambda x
: self
.from_vector(x
), ls
))
118 def _element_constructor_(self
, elt
):
120 Construct an element of this algebra from its natural
123 This gets called only after the parent element _call_ method
124 fails to find a coercion for the argument.
128 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
130 ....: RealSymmetricEJA)
134 The identity in `S^n` is converted to the identity in the EJA::
136 sage: J = RealSymmetricEJA(3)
137 sage: I = matrix.identity(QQ,3)
138 sage: J(I) == J.one()
141 This skew-symmetric matrix can't be represented in the EJA::
143 sage: J = RealSymmetricEJA(3)
144 sage: A = matrix(QQ,3, lambda i,j: i-j)
146 Traceback (most recent call last):
148 ArithmeticError: vector is not in free module
152 Ensure that we can convert any element of the two non-matrix
153 simple algebras (whose natural representations are their usual
154 vector representations) back and forth faithfully::
156 sage: set_random_seed()
157 sage: J = HadamardEJA.random_instance()
158 sage: x = J.random_element()
159 sage: J(x.to_vector().column()) == x
161 sage: J = JordanSpinEJA.random_instance()
162 sage: x = J.random_element()
163 sage: J(x.to_vector().column()) == x
167 msg
= "not a naturally-represented algebra element"
169 # The superclass implementation of random_element()
170 # needs to be able to coerce "0" into the algebra.
172 elif elt
in self
.base_ring():
173 # Ensure that no base ring -> algebra coercion is performed
174 # by this method. There's some stupidity in sage that would
175 # otherwise propagate to this method; for example, sage thinks
176 # that the integer 3 belongs to the space of 2-by-2 matrices.
177 raise ValueError(msg
)
179 natural_basis
= self
.natural_basis()
180 basis_space
= natural_basis
[0].matrix_space()
181 if elt
not in basis_space
:
182 raise ValueError(msg
)
184 # Thanks for nothing! Matrix spaces aren't vector spaces in
185 # Sage, so we have to figure out its natural-basis coordinates
186 # ourselves. We use the basis space's ring instead of the
187 # element's ring because the basis space might be an algebraic
188 # closure whereas the base ring of the 3-by-3 identity matrix
189 # could be QQ instead of QQbar.
190 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
191 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
192 coords
= W
.coordinate_vector(_mat2vec(elt
))
193 return self
.from_vector(coords
)
196 def _max_test_case_size():
198 Return an integer "size" that is an upper bound on the size of
199 this algebra when it is used in a random test
200 case. Unfortunately, the term "size" is quite vague -- when
201 dealing with `R^n` under either the Hadamard or Jordan spin
202 product, the "size" refers to the dimension `n`. When dealing
203 with a matrix algebra (real symmetric or complex/quaternion
204 Hermitian), it refers to the size of the matrix, which is
205 far less than the dimension of the underlying vector space.
207 We default to five in this class, which is safe in `R^n`. The
208 matrix algebra subclasses (or any class where the "size" is
209 interpreted to be far less than the dimension) should override
210 with a smaller number.
216 Return a string representation of ``self``.
220 sage: from mjo.eja.eja_algebra import JordanSpinEJA
224 Ensure that it says what we think it says::
226 sage: JordanSpinEJA(2, field=AA)
227 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
228 sage: JordanSpinEJA(3, field=RDF)
229 Euclidean Jordan algebra of dimension 3 over Real Double Field
232 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
233 return fmt
.format(self
.dimension(), self
.base_ring())
235 def product_on_basis(self
, i
, j
):
236 return self
._multiplication
_table
[i
][j
]
239 def characteristic_polynomial_of(self
):
241 Return the algebra's "characteristic polynomial of" function,
242 which is itself a multivariate polynomial that, when evaluated
243 at the coordinates of some algebra element, returns that
244 element's characteristic polynomial.
246 The resulting polynomial has `n+1` variables, where `n` is the
247 dimension of this algebra. The first `n` variables correspond to
248 the coordinates of an algebra element: when evaluated at the
249 coordinates of an algebra element with respect to a certain
250 basis, the result is a univariate polynomial (in the one
251 remaining variable ``t``), namely the characteristic polynomial
256 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
260 The characteristic polynomial in the spin algebra is given in
261 Alizadeh, Example 11.11::
263 sage: J = JordanSpinEJA(3)
264 sage: p = J.characteristic_polynomial_of(); p
265 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
266 sage: xvec = J.one().to_vector()
270 By definition, the characteristic polynomial is a monic
271 degree-zero polynomial in a rank-zero algebra. Note that
272 Cayley-Hamilton is indeed satisfied since the polynomial
273 ``1`` evaluates to the identity element of the algebra on
276 sage: J = TrivialEJA()
277 sage: J.characteristic_polynomial_of()
284 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
285 a
= self
._charpoly
_coefficients
()
287 # We go to a bit of trouble here to reorder the
288 # indeterminates, so that it's easier to evaluate the
289 # characteristic polynomial at x's coordinates and get back
290 # something in terms of t, which is what we want.
291 S
= PolynomialRing(self
.base_ring(),'t')
295 S
= PolynomialRing(S
, R
.variable_names())
298 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
301 def inner_product(self
, x
, y
):
303 The inner product associated with this Euclidean Jordan algebra.
305 Defaults to the trace inner product, but can be overridden by
306 subclasses if they are sure that the necessary properties are
311 sage: from mjo.eja.eja_algebra import random_eja
315 Our inner product is "associative," which means the following for
316 a symmetric bilinear form::
318 sage: set_random_seed()
319 sage: J = random_eja()
320 sage: x,y,z = J.random_elements(3)
321 sage: (x*y).inner_product(z) == y.inner_product(x*z)
325 X
= x
.natural_representation()
326 Y
= y
.natural_representation()
327 return self
.natural_inner_product(X
,Y
)
330 def is_trivial(self
):
332 Return whether or not this algebra is trivial.
334 A trivial algebra contains only the zero element.
338 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
343 sage: J = ComplexHermitianEJA(3)
349 sage: J = TrivialEJA()
354 return self
.dimension() == 0
357 def multiplication_table(self
):
359 Return a visual representation of this algebra's multiplication
360 table (on basis elements).
364 sage: from mjo.eja.eja_algebra import JordanSpinEJA
368 sage: J = JordanSpinEJA(4)
369 sage: J.multiplication_table()
370 +----++----+----+----+----+
371 | * || e0 | e1 | e2 | e3 |
372 +====++====+====+====+====+
373 | e0 || e0 | e1 | e2 | e3 |
374 +----++----+----+----+----+
375 | e1 || e1 | e0 | 0 | 0 |
376 +----++----+----+----+----+
377 | e2 || e2 | 0 | e0 | 0 |
378 +----++----+----+----+----+
379 | e3 || e3 | 0 | 0 | e0 |
380 +----++----+----+----+----+
383 M
= list(self
._multiplication
_table
) # copy
384 for i
in range(len(M
)):
385 # M had better be "square"
386 M
[i
] = [self
.monomial(i
)] + M
[i
]
387 M
= [["*"] + list(self
.gens())] + M
388 return table(M
, header_row
=True, header_column
=True, frame
=True)
391 def natural_basis(self
):
393 Return a more-natural representation of this algebra's basis.
395 Every finite-dimensional Euclidean Jordan Algebra is a direct
396 sum of five simple algebras, four of which comprise Hermitian
397 matrices. This method returns the original "natural" basis
398 for our underlying vector space. (Typically, the natural basis
399 is used to construct the multiplication table in the first place.)
401 Note that this will always return a matrix. The standard basis
402 in `R^n` will be returned as `n`-by-`1` column matrices.
406 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
407 ....: RealSymmetricEJA)
411 sage: J = RealSymmetricEJA(2)
413 Finite family {0: e0, 1: e1, 2: e2}
414 sage: J.natural_basis()
416 [1 0] [ 0 0.7071067811865475?] [0 0]
417 [0 0], [0.7071067811865475? 0], [0 1]
422 sage: J = JordanSpinEJA(2)
424 Finite family {0: e0, 1: e1}
425 sage: J.natural_basis()
432 if self
._natural
_basis
is None:
433 M
= self
.natural_basis_space()
434 return tuple( M(b
.to_vector()) for b
in self
.basis() )
436 return self
._natural
_basis
439 def natural_basis_space(self
):
441 Return the matrix space in which this algebra's natural basis
444 Generally this will be an `n`-by-`1` column-vector space,
445 except when the algebra is trivial. There it's `n`-by-`n`
446 (where `n` is zero), to ensure that two elements of the
447 natural basis space (empty matrices) can be multiplied.
449 if self
.is_trivial():
450 return MatrixSpace(self
.base_ring(), 0)
451 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
452 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
454 return self
._natural
_basis
[0].matrix_space()
458 def natural_inner_product(X
,Y
):
460 Compute the inner product of two naturally-represented elements.
462 For example in the real symmetric matrix EJA, this will compute
463 the trace inner-product of two n-by-n symmetric matrices. The
464 default should work for the real cartesian product EJA, the
465 Jordan spin EJA, and the real symmetric matrices. The others
466 will have to be overridden.
468 return (X
.conjugate_transpose()*Y
).trace()
474 Return the unit element of this algebra.
478 sage: from mjo.eja.eja_algebra import (HadamardEJA,
483 sage: J = HadamardEJA(5)
485 e0 + e1 + e2 + e3 + e4
489 The identity element acts like the identity::
491 sage: set_random_seed()
492 sage: J = random_eja()
493 sage: x = J.random_element()
494 sage: J.one()*x == x and x*J.one() == x
497 The matrix of the unit element's operator is the identity::
499 sage: set_random_seed()
500 sage: J = random_eja()
501 sage: actual = J.one().operator().matrix()
502 sage: expected = matrix.identity(J.base_ring(), J.dimension())
503 sage: actual == expected
507 # We can brute-force compute the matrices of the operators
508 # that correspond to the basis elements of this algebra.
509 # If some linear combination of those basis elements is the
510 # algebra identity, then the same linear combination of
511 # their matrices has to be the identity matrix.
513 # Of course, matrices aren't vectors in sage, so we have to
514 # appeal to the "long vectors" isometry.
515 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
517 # Now we use basis linear algebra to find the coefficients,
518 # of the matrices-as-vectors-linear-combination, which should
519 # work for the original algebra basis too.
520 A
= matrix
.column(self
.base_ring(), oper_vecs
)
522 # We used the isometry on the left-hand side already, but we
523 # still need to do it for the right-hand side. Recall that we
524 # wanted something that summed to the identity matrix.
525 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
527 # Now if there's an identity element in the algebra, this should work.
528 coeffs
= A
.solve_right(b
)
529 return self
.linear_combination(zip(self
.gens(), coeffs
))
532 def peirce_decomposition(self
, c
):
534 The Peirce decomposition of this algebra relative to the
537 In the future, this can be extended to a complete system of
538 orthogonal idempotents.
542 - ``c`` -- an idempotent of this algebra.
546 A triple (J0, J5, J1) containing two subalgebras and one subspace
549 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
550 corresponding to the eigenvalue zero.
552 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
553 corresponding to the eigenvalue one-half.
555 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
556 corresponding to the eigenvalue one.
558 These are the only possible eigenspaces for that operator, and this
559 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
560 orthogonal, and are subalgebras of this algebra with the appropriate
565 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
569 The canonical example comes from the symmetric matrices, which
570 decompose into diagonal and off-diagonal parts::
572 sage: J = RealSymmetricEJA(3)
573 sage: C = matrix(QQ, [ [1,0,0],
577 sage: J0,J5,J1 = J.peirce_decomposition(c)
579 Euclidean Jordan algebra of dimension 1...
581 Vector space of degree 6 and dimension 2...
583 Euclidean Jordan algebra of dimension 3...
584 sage: J0.one().natural_representation()
588 sage: orig_df = AA.options.display_format
589 sage: AA.options.display_format = 'radical'
590 sage: J.from_vector(J5.basis()[0]).natural_representation()
594 sage: J.from_vector(J5.basis()[1]).natural_representation()
598 sage: AA.options.display_format = orig_df
599 sage: J1.one().natural_representation()
606 Every algebra decomposes trivially with respect to its identity
609 sage: set_random_seed()
610 sage: J = random_eja()
611 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
612 sage: J0.dimension() == 0 and J5.dimension() == 0
614 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
617 The decomposition is into eigenspaces, and its components are
618 therefore necessarily orthogonal. Moreover, the identity
619 elements in the two subalgebras are the projections onto their
620 respective subspaces of the superalgebra's identity element::
622 sage: set_random_seed()
623 sage: J = random_eja()
624 sage: x = J.random_element()
625 sage: if not J.is_trivial():
626 ....: while x.is_nilpotent():
627 ....: x = J.random_element()
628 sage: c = x.subalgebra_idempotent()
629 sage: J0,J5,J1 = J.peirce_decomposition(c)
631 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
632 ....: w = w.superalgebra_element()
633 ....: y = J.from_vector(y)
634 ....: z = z.superalgebra_element()
635 ....: ipsum += w.inner_product(y).abs()
636 ....: ipsum += w.inner_product(z).abs()
637 ....: ipsum += y.inner_product(z).abs()
640 sage: J1(c) == J1.one()
642 sage: J0(J.one() - c) == J0.one()
646 if not c
.is_idempotent():
647 raise ValueError("element is not idempotent: %s" % c
)
649 # Default these to what they should be if they turn out to be
650 # trivial, because eigenspaces_left() won't return eigenvalues
651 # corresponding to trivial spaces (e.g. it returns only the
652 # eigenspace corresponding to lambda=1 if you take the
653 # decomposition relative to the identity element).
654 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
655 J0
= trivial
# eigenvalue zero
656 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
657 J1
= trivial
# eigenvalue one
659 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
660 if eigval
== ~
(self
.base_ring()(2)):
663 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
664 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
670 raise ValueError("unexpected eigenvalue: %s" % eigval
)
675 def random_elements(self
, count
):
677 Return ``count`` random elements as a tuple.
681 sage: from mjo.eja.eja_algebra import JordanSpinEJA
685 sage: J = JordanSpinEJA(3)
686 sage: x,y,z = J.random_elements(3)
687 sage: all( [ x in J, y in J, z in J ])
689 sage: len( J.random_elements(10) ) == 10
693 return tuple( self
.random_element() for idx
in range(count
) )
696 def random_instance(cls
, field
=AA
, **kwargs
):
698 Return a random instance of this type of algebra.
700 Beware, this will crash for "most instances" because the
701 constructor below looks wrong.
703 if cls
is TrivialEJA
:
704 # The TrivialEJA class doesn't take an "n" argument because
708 n
= ZZ
.random_element(cls
._max
_test
_case
_size
() + 1)
709 return cls(n
, field
, **kwargs
)
712 def _charpoly_coefficients(self
):
714 The `r` polynomial coefficients of the "characteristic polynomial
718 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
719 R
= PolynomialRing(self
.base_ring(), var_names
)
721 F
= R
.fraction_field()
724 # From a result in my book, these are the entries of the
725 # basis representation of L_x.
726 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
729 L_x
= matrix(F
, n
, n
, L_x_i_j
)
732 if self
.rank
.is_in_cache():
734 # There's no need to pad the system with redundant
735 # columns if we *know* they'll be redundant.
738 # Compute an extra power in case the rank is equal to
739 # the dimension (otherwise, we would stop at x^(r-1)).
740 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
741 for k
in range(n
+1) ]
742 A
= matrix
.column(F
, x_powers
[:n
])
743 AE
= A
.extended_echelon_form()
750 # The theory says that only the first "r" coefficients are
751 # nonzero, and they actually live in the original polynomial
752 # ring and not the fraction field. We negate them because
753 # in the actual characteristic polynomial, they get moved
754 # to the other side where x^r lives.
755 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
760 Return the rank of this EJA.
762 This is a cached method because we know the rank a priori for
763 all of the algebras we can construct. Thus we can avoid the
764 expensive ``_charpoly_coefficients()`` call unless we truly
765 need to compute the whole characteristic polynomial.
769 sage: from mjo.eja.eja_algebra import (HadamardEJA,
771 ....: RealSymmetricEJA,
772 ....: ComplexHermitianEJA,
773 ....: QuaternionHermitianEJA,
778 The rank of the Jordan spin algebra is always two::
780 sage: JordanSpinEJA(2).rank()
782 sage: JordanSpinEJA(3).rank()
784 sage: JordanSpinEJA(4).rank()
787 The rank of the `n`-by-`n` Hermitian real, complex, or
788 quaternion matrices is `n`::
790 sage: RealSymmetricEJA(4).rank()
792 sage: ComplexHermitianEJA(3).rank()
794 sage: QuaternionHermitianEJA(2).rank()
799 Ensure that every EJA that we know how to construct has a
800 positive integer rank, unless the algebra is trivial in
801 which case its rank will be zero::
803 sage: set_random_seed()
804 sage: J = random_eja()
808 sage: r > 0 or (r == 0 and J.is_trivial())
811 Ensure that computing the rank actually works, since the ranks
812 of all simple algebras are known and will be cached by default::
814 sage: J = HadamardEJA(4)
815 sage: J.rank.clear_cache()
821 sage: J = JordanSpinEJA(4)
822 sage: J.rank.clear_cache()
828 sage: J = RealSymmetricEJA(3)
829 sage: J.rank.clear_cache()
835 sage: J = ComplexHermitianEJA(2)
836 sage: J.rank.clear_cache()
842 sage: J = QuaternionHermitianEJA(2)
843 sage: J.rank.clear_cache()
847 return len(self
._charpoly
_coefficients
())
850 def vector_space(self
):
852 Return the vector space that underlies this algebra.
856 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
860 sage: J = RealSymmetricEJA(2)
861 sage: J.vector_space()
862 Vector space of dimension 3 over...
865 return self
.zero().to_vector().parent().ambient_vector_space()
868 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
871 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
):
873 Return the Euclidean Jordan Algebra corresponding to the set
874 `R^n` under the Hadamard product.
876 Note: this is nothing more than the Cartesian product of ``n``
877 copies of the spin algebra. Once Cartesian product algebras
878 are implemented, this can go.
882 sage: from mjo.eja.eja_algebra import HadamardEJA
886 This multiplication table can be verified by hand::
888 sage: J = HadamardEJA(3)
889 sage: e0,e1,e2 = J.gens()
905 We can change the generator prefix::
907 sage: HadamardEJA(3, prefix='r').gens()
911 def __init__(self
, n
, field
=AA
, **kwargs
):
912 V
= VectorSpace(field
, n
)
913 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
916 fdeja
= super(HadamardEJA
, self
)
917 fdeja
.__init
__(field
, mult_table
, **kwargs
)
918 self
.rank
.set_cache(n
)
920 def inner_product(self
, x
, y
):
922 Faster to reimplement than to use natural representations.
926 sage: from mjo.eja.eja_algebra import HadamardEJA
930 Ensure that this is the usual inner product for the algebras
933 sage: set_random_seed()
934 sage: J = HadamardEJA.random_instance()
935 sage: x,y = J.random_elements(2)
936 sage: X = x.natural_representation()
937 sage: Y = y.natural_representation()
938 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
942 return x
.to_vector().inner_product(y
.to_vector())
945 def random_eja(field
=AA
):
947 Return a "random" finite-dimensional Euclidean Jordan Algebra.
951 sage: from mjo.eja.eja_algebra import random_eja
956 Euclidean Jordan algebra of dimension...
959 classname
= choice([TrivialEJA
,
964 QuaternionHermitianEJA
])
965 return classname
.random_instance(field
=field
)
970 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
972 def _max_test_case_size():
973 # Play it safe, since this will be squared and the underlying
974 # field can have dimension 4 (quaternions) too.
977 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
979 Compared to the superclass constructor, we take a basis instead of
980 a multiplication table because the latter can be computed in terms
981 of the former when the product is known (like it is here).
983 # Used in this class's fast _charpoly_coefficients() override.
984 self
._basis
_normalizers
= None
986 # We're going to loop through this a few times, so now's a good
987 # time to ensure that it isn't a generator expression.
990 if len(basis
) > 1 and normalize_basis
:
991 # We'll need sqrt(2) to normalize the basis, and this
992 # winds up in the multiplication table, so the whole
993 # algebra needs to be over the field extension.
994 R
= PolynomialRing(field
, 'z')
997 if p
.is_irreducible():
998 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
999 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1000 self
._basis
_normalizers
= tuple(
1001 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1002 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1004 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1006 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
1007 fdeja
.__init
__(field
, Qs
, natural_basis
=basis
, **kwargs
)
1012 def _charpoly_coefficients(self
):
1014 Override the parent method with something that tries to compute
1015 over a faster (non-extension) field.
1017 if self
._basis
_normalizers
is None:
1018 # We didn't normalize, so assume that the basis we started
1019 # with had entries in a nice field.
1020 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1022 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1023 self
._basis
_normalizers
) )
1025 # Do this over the rationals and convert back at the end.
1026 # Only works because we know the entries of the basis are
1028 J
= MatrixEuclideanJordanAlgebra(QQ
,
1030 normalize_basis
=False)
1031 a
= J
._charpoly
_coefficients
()
1033 # Unfortunately, changing the basis does change the
1034 # coefficients of the characteristic polynomial, but since
1035 # these are really the coefficients of the "characteristic
1036 # polynomial of" function, everything is still nice and
1037 # unevaluated. It's therefore "obvious" how scaling the
1038 # basis affects the coordinate variables X1, X2, et
1039 # cetera. Scaling the first basis vector up by "n" adds a
1040 # factor of 1/n into every "X1" term, for example. So here
1041 # we simply undo the basis_normalizer scaling that we
1042 # performed earlier.
1044 # The a[0] access here is safe because trivial algebras
1045 # won't have any basis normalizers and therefore won't
1046 # make it to this "else" branch.
1047 XS
= a
[0].parent().gens()
1048 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1049 for i
in range(len(XS
)) }
1050 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1054 def multiplication_table_from_matrix_basis(basis
):
1056 At least three of the five simple Euclidean Jordan algebras have the
1057 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1058 multiplication on the right is matrix multiplication. Given a basis
1059 for the underlying matrix space, this function returns a
1060 multiplication table (obtained by looping through the basis
1061 elements) for an algebra of those matrices.
1063 # In S^2, for example, we nominally have four coordinates even
1064 # though the space is of dimension three only. The vector space V
1065 # is supposed to hold the entire long vector, and the subspace W
1066 # of V will be spanned by the vectors that arise from symmetric
1067 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1071 field
= basis
[0].base_ring()
1072 dimension
= basis
[0].nrows()
1074 V
= VectorSpace(field
, dimension
**2)
1075 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1077 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1080 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1081 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1089 Embed the matrix ``M`` into a space of real matrices.
1091 The matrix ``M`` can have entries in any field at the moment:
1092 the real numbers, complex numbers, or quaternions. And although
1093 they are not a field, we can probably support octonions at some
1094 point, too. This function returns a real matrix that "acts like"
1095 the original with respect to matrix multiplication; i.e.
1097 real_embed(M*N) = real_embed(M)*real_embed(N)
1100 raise NotImplementedError
1104 def real_unembed(M
):
1106 The inverse of :meth:`real_embed`.
1108 raise NotImplementedError
1112 def natural_inner_product(cls
,X
,Y
):
1113 Xu
= cls
.real_unembed(X
)
1114 Yu
= cls
.real_unembed(Y
)
1115 tr
= (Xu
*Yu
).trace()
1118 # It's real already.
1121 # Otherwise, try the thing that works for complex numbers; and
1122 # if that doesn't work, the thing that works for quaternions.
1124 return tr
.vector()[0] # real part, imag part is index 1
1125 except AttributeError:
1126 # A quaternions doesn't have a vector() method, but does
1127 # have coefficient_tuple() method that returns the
1128 # coefficients of 1, i, j, and k -- in that order.
1129 return tr
.coefficient_tuple()[0]
1132 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1136 The identity function, for embedding real matrices into real
1142 def real_unembed(M
):
1144 The identity function, for unembedding real matrices from real
1150 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1152 The rank-n simple EJA consisting of real symmetric n-by-n
1153 matrices, the usual symmetric Jordan product, and the trace inner
1154 product. It has dimension `(n^2 + n)/2` over the reals.
1158 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1162 sage: J = RealSymmetricEJA(2)
1163 sage: e0, e1, e2 = J.gens()
1171 In theory, our "field" can be any subfield of the reals::
1173 sage: RealSymmetricEJA(2, RDF)
1174 Euclidean Jordan algebra of dimension 3 over Real Double Field
1175 sage: RealSymmetricEJA(2, RR)
1176 Euclidean Jordan algebra of dimension 3 over Real Field with
1177 53 bits of precision
1181 The dimension of this algebra is `(n^2 + n) / 2`::
1183 sage: set_random_seed()
1184 sage: n_max = RealSymmetricEJA._max_test_case_size()
1185 sage: n = ZZ.random_element(1, n_max)
1186 sage: J = RealSymmetricEJA(n)
1187 sage: J.dimension() == (n^2 + n)/2
1190 The Jordan multiplication is what we think it is::
1192 sage: set_random_seed()
1193 sage: J = RealSymmetricEJA.random_instance()
1194 sage: x,y = J.random_elements(2)
1195 sage: actual = (x*y).natural_representation()
1196 sage: X = x.natural_representation()
1197 sage: Y = y.natural_representation()
1198 sage: expected = (X*Y + Y*X)/2
1199 sage: actual == expected
1201 sage: J(expected) == x*y
1204 We can change the generator prefix::
1206 sage: RealSymmetricEJA(3, prefix='q').gens()
1207 (q0, q1, q2, q3, q4, q5)
1209 Our natural basis is normalized with respect to the natural inner
1210 product unless we specify otherwise::
1212 sage: set_random_seed()
1213 sage: J = RealSymmetricEJA.random_instance()
1214 sage: all( b.norm() == 1 for b in J.gens() )
1217 Since our natural basis is normalized with respect to the natural
1218 inner product, and since we know that this algebra is an EJA, any
1219 left-multiplication operator's matrix will be symmetric because
1220 natural->EJA basis representation is an isometry and within the EJA
1221 the operator is self-adjoint by the Jordan axiom::
1223 sage: set_random_seed()
1224 sage: x = RealSymmetricEJA.random_instance().random_element()
1225 sage: x.operator().matrix().is_symmetric()
1228 We can construct the (trivial) algebra of rank zero::
1230 sage: RealSymmetricEJA(0)
1231 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1235 def _denormalized_basis(cls
, n
, field
):
1237 Return a basis for the space of real symmetric n-by-n matrices.
1241 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1245 sage: set_random_seed()
1246 sage: n = ZZ.random_element(1,5)
1247 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1248 sage: all( M.is_symmetric() for M in B)
1252 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1256 for j
in range(i
+1):
1257 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1261 Sij
= Eij
+ Eij
.transpose()
1267 def _max_test_case_size():
1268 return 4 # Dimension 10
1271 def __init__(self
, n
, field
=AA
, **kwargs
):
1272 basis
= self
._denormalized
_basis
(n
, field
)
1273 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, **kwargs
)
1274 self
.rank
.set_cache(n
)
1277 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1281 Embed the n-by-n complex matrix ``M`` into the space of real
1282 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1283 bi` to the block matrix ``[[a,b],[-b,a]]``.
1287 sage: from mjo.eja.eja_algebra import \
1288 ....: ComplexMatrixEuclideanJordanAlgebra
1292 sage: F = QuadraticField(-1, 'I')
1293 sage: x1 = F(4 - 2*i)
1294 sage: x2 = F(1 + 2*i)
1297 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1298 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1307 Embedding is a homomorphism (isomorphism, in fact)::
1309 sage: set_random_seed()
1310 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1311 sage: n = ZZ.random_element(n_max)
1312 sage: F = QuadraticField(-1, 'I')
1313 sage: X = random_matrix(F, n)
1314 sage: Y = random_matrix(F, n)
1315 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1316 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1317 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1324 raise ValueError("the matrix 'M' must be square")
1326 # We don't need any adjoined elements...
1327 field
= M
.base_ring().base_ring()
1331 a
= z
.list()[0] # real part, I guess
1332 b
= z
.list()[1] # imag part, I guess
1333 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1335 return matrix
.block(field
, n
, blocks
)
1339 def real_unembed(M
):
1341 The inverse of _embed_complex_matrix().
1345 sage: from mjo.eja.eja_algebra import \
1346 ....: ComplexMatrixEuclideanJordanAlgebra
1350 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1351 ....: [-2, 1, -4, 3],
1352 ....: [ 9, 10, 11, 12],
1353 ....: [-10, 9, -12, 11] ])
1354 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1356 [ 10*I + 9 12*I + 11]
1360 Unembedding is the inverse of embedding::
1362 sage: set_random_seed()
1363 sage: F = QuadraticField(-1, 'I')
1364 sage: M = random_matrix(F, 3)
1365 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1366 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1372 raise ValueError("the matrix 'M' must be square")
1373 if not n
.mod(2).is_zero():
1374 raise ValueError("the matrix 'M' must be a complex embedding")
1376 # If "M" was normalized, its base ring might have roots
1377 # adjoined and they can stick around after unembedding.
1378 field
= M
.base_ring()
1379 R
= PolynomialRing(field
, 'z')
1382 # Sage doesn't know how to embed AA into QQbar, i.e. how
1383 # to adjoin sqrt(-1) to AA.
1386 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1389 # Go top-left to bottom-right (reading order), converting every
1390 # 2-by-2 block we see to a single complex element.
1392 for k
in range(n
/2):
1393 for j
in range(n
/2):
1394 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1395 if submat
[0,0] != submat
[1,1]:
1396 raise ValueError('bad on-diagonal submatrix')
1397 if submat
[0,1] != -submat
[1,0]:
1398 raise ValueError('bad off-diagonal submatrix')
1399 z
= submat
[0,0] + submat
[0,1]*i
1402 return matrix(F
, n
/2, elements
)
1406 def natural_inner_product(cls
,X
,Y
):
1408 Compute a natural inner product in this algebra directly from
1413 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1417 This gives the same answer as the slow, default method implemented
1418 in :class:`MatrixEuclideanJordanAlgebra`::
1420 sage: set_random_seed()
1421 sage: J = ComplexHermitianEJA.random_instance()
1422 sage: x,y = J.random_elements(2)
1423 sage: Xe = x.natural_representation()
1424 sage: Ye = y.natural_representation()
1425 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1426 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1427 sage: expected = (X*Y).trace().real()
1428 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1429 sage: actual == expected
1433 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1436 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1438 The rank-n simple EJA consisting of complex Hermitian n-by-n
1439 matrices over the real numbers, the usual symmetric Jordan product,
1440 and the real-part-of-trace inner product. It has dimension `n^2` over
1445 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1449 In theory, our "field" can be any subfield of the reals::
1451 sage: ComplexHermitianEJA(2, RDF)
1452 Euclidean Jordan algebra of dimension 4 over Real Double Field
1453 sage: ComplexHermitianEJA(2, RR)
1454 Euclidean Jordan algebra of dimension 4 over Real Field with
1455 53 bits of precision
1459 The dimension of this algebra is `n^2`::
1461 sage: set_random_seed()
1462 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1463 sage: n = ZZ.random_element(1, n_max)
1464 sage: J = ComplexHermitianEJA(n)
1465 sage: J.dimension() == n^2
1468 The Jordan multiplication is what we think it is::
1470 sage: set_random_seed()
1471 sage: J = ComplexHermitianEJA.random_instance()
1472 sage: x,y = J.random_elements(2)
1473 sage: actual = (x*y).natural_representation()
1474 sage: X = x.natural_representation()
1475 sage: Y = y.natural_representation()
1476 sage: expected = (X*Y + Y*X)/2
1477 sage: actual == expected
1479 sage: J(expected) == x*y
1482 We can change the generator prefix::
1484 sage: ComplexHermitianEJA(2, prefix='z').gens()
1487 Our natural basis is normalized with respect to the natural inner
1488 product unless we specify otherwise::
1490 sage: set_random_seed()
1491 sage: J = ComplexHermitianEJA.random_instance()
1492 sage: all( b.norm() == 1 for b in J.gens() )
1495 Since our natural basis is normalized with respect to the natural
1496 inner product, and since we know that this algebra is an EJA, any
1497 left-multiplication operator's matrix will be symmetric because
1498 natural->EJA basis representation is an isometry and within the EJA
1499 the operator is self-adjoint by the Jordan axiom::
1501 sage: set_random_seed()
1502 sage: x = ComplexHermitianEJA.random_instance().random_element()
1503 sage: x.operator().matrix().is_symmetric()
1506 We can construct the (trivial) algebra of rank zero::
1508 sage: ComplexHermitianEJA(0)
1509 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1514 def _denormalized_basis(cls
, n
, field
):
1516 Returns a basis for the space of complex Hermitian n-by-n matrices.
1518 Why do we embed these? Basically, because all of numerical linear
1519 algebra assumes that you're working with vectors consisting of `n`
1520 entries from a field and scalars from the same field. There's no way
1521 to tell SageMath that (for example) the vectors contain complex
1522 numbers, while the scalar field is real.
1526 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1530 sage: set_random_seed()
1531 sage: n = ZZ.random_element(1,5)
1532 sage: field = QuadraticField(2, 'sqrt2')
1533 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1534 sage: all( M.is_symmetric() for M in B)
1538 R
= PolynomialRing(field
, 'z')
1540 F
= field
.extension(z
**2 + 1, 'I')
1543 # This is like the symmetric case, but we need to be careful:
1545 # * We want conjugate-symmetry, not just symmetry.
1546 # * The diagonal will (as a result) be real.
1550 for j
in range(i
+1):
1551 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1553 Sij
= cls
.real_embed(Eij
)
1556 # The second one has a minus because it's conjugated.
1557 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1559 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1562 # Since we embedded these, we can drop back to the "field" that we
1563 # started with instead of the complex extension "F".
1564 return ( s
.change_ring(field
) for s
in S
)
1567 def __init__(self
, n
, field
=AA
, **kwargs
):
1568 basis
= self
._denormalized
_basis
(n
,field
)
1569 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1570 self
.rank
.set_cache(n
)
1573 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1577 Embed the n-by-n quaternion matrix ``M`` into the space of real
1578 matrices of size 4n-by-4n by first sending each quaternion entry `z
1579 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1580 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1585 sage: from mjo.eja.eja_algebra import \
1586 ....: QuaternionMatrixEuclideanJordanAlgebra
1590 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1591 sage: i,j,k = Q.gens()
1592 sage: x = 1 + 2*i + 3*j + 4*k
1593 sage: M = matrix(Q, 1, [[x]])
1594 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1600 Embedding is a homomorphism (isomorphism, in fact)::
1602 sage: set_random_seed()
1603 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1604 sage: n = ZZ.random_element(n_max)
1605 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1606 sage: X = random_matrix(Q, n)
1607 sage: Y = random_matrix(Q, n)
1608 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1609 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1610 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1615 quaternions
= M
.base_ring()
1618 raise ValueError("the matrix 'M' must be square")
1620 F
= QuadraticField(-1, 'I')
1625 t
= z
.coefficient_tuple()
1630 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1631 [-c
+ d
*i
, a
- b
*i
]])
1632 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1633 blocks
.append(realM
)
1635 # We should have real entries by now, so use the realest field
1636 # we've got for the return value.
1637 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1642 def real_unembed(M
):
1644 The inverse of _embed_quaternion_matrix().
1648 sage: from mjo.eja.eja_algebra import \
1649 ....: QuaternionMatrixEuclideanJordanAlgebra
1653 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1654 ....: [-2, 1, -4, 3],
1655 ....: [-3, 4, 1, -2],
1656 ....: [-4, -3, 2, 1]])
1657 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1658 [1 + 2*i + 3*j + 4*k]
1662 Unembedding is the inverse of embedding::
1664 sage: set_random_seed()
1665 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1666 sage: M = random_matrix(Q, 3)
1667 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1668 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1674 raise ValueError("the matrix 'M' must be square")
1675 if not n
.mod(4).is_zero():
1676 raise ValueError("the matrix 'M' must be a quaternion embedding")
1678 # Use the base ring of the matrix to ensure that its entries can be
1679 # multiplied by elements of the quaternion algebra.
1680 field
= M
.base_ring()
1681 Q
= QuaternionAlgebra(field
,-1,-1)
1684 # Go top-left to bottom-right (reading order), converting every
1685 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1688 for l
in range(n
/4):
1689 for m
in range(n
/4):
1690 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1691 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1692 if submat
[0,0] != submat
[1,1].conjugate():
1693 raise ValueError('bad on-diagonal submatrix')
1694 if submat
[0,1] != -submat
[1,0].conjugate():
1695 raise ValueError('bad off-diagonal submatrix')
1696 z
= submat
[0,0].real()
1697 z
+= submat
[0,0].imag()*i
1698 z
+= submat
[0,1].real()*j
1699 z
+= submat
[0,1].imag()*k
1702 return matrix(Q
, n
/4, elements
)
1706 def natural_inner_product(cls
,X
,Y
):
1708 Compute a natural inner product in this algebra directly from
1713 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1717 This gives the same answer as the slow, default method implemented
1718 in :class:`MatrixEuclideanJordanAlgebra`::
1720 sage: set_random_seed()
1721 sage: J = QuaternionHermitianEJA.random_instance()
1722 sage: x,y = J.random_elements(2)
1723 sage: Xe = x.natural_representation()
1724 sage: Ye = y.natural_representation()
1725 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1726 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1727 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1728 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1729 sage: actual == expected
1733 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1736 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1738 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1739 matrices, the usual symmetric Jordan product, and the
1740 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1745 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1749 In theory, our "field" can be any subfield of the reals::
1751 sage: QuaternionHermitianEJA(2, RDF)
1752 Euclidean Jordan algebra of dimension 6 over Real Double Field
1753 sage: QuaternionHermitianEJA(2, RR)
1754 Euclidean Jordan algebra of dimension 6 over Real Field with
1755 53 bits of precision
1759 The dimension of this algebra is `2*n^2 - n`::
1761 sage: set_random_seed()
1762 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1763 sage: n = ZZ.random_element(1, n_max)
1764 sage: J = QuaternionHermitianEJA(n)
1765 sage: J.dimension() == 2*(n^2) - n
1768 The Jordan multiplication is what we think it is::
1770 sage: set_random_seed()
1771 sage: J = QuaternionHermitianEJA.random_instance()
1772 sage: x,y = J.random_elements(2)
1773 sage: actual = (x*y).natural_representation()
1774 sage: X = x.natural_representation()
1775 sage: Y = y.natural_representation()
1776 sage: expected = (X*Y + Y*X)/2
1777 sage: actual == expected
1779 sage: J(expected) == x*y
1782 We can change the generator prefix::
1784 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1785 (a0, a1, a2, a3, a4, a5)
1787 Our natural basis is normalized with respect to the natural inner
1788 product unless we specify otherwise::
1790 sage: set_random_seed()
1791 sage: J = QuaternionHermitianEJA.random_instance()
1792 sage: all( b.norm() == 1 for b in J.gens() )
1795 Since our natural basis is normalized with respect to the natural
1796 inner product, and since we know that this algebra is an EJA, any
1797 left-multiplication operator's matrix will be symmetric because
1798 natural->EJA basis representation is an isometry and within the EJA
1799 the operator is self-adjoint by the Jordan axiom::
1801 sage: set_random_seed()
1802 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1803 sage: x.operator().matrix().is_symmetric()
1806 We can construct the (trivial) algebra of rank zero::
1808 sage: QuaternionHermitianEJA(0)
1809 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1813 def _denormalized_basis(cls
, n
, field
):
1815 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1817 Why do we embed these? Basically, because all of numerical
1818 linear algebra assumes that you're working with vectors consisting
1819 of `n` entries from a field and scalars from the same field. There's
1820 no way to tell SageMath that (for example) the vectors contain
1821 complex numbers, while the scalar field is real.
1825 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1829 sage: set_random_seed()
1830 sage: n = ZZ.random_element(1,5)
1831 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1832 sage: all( M.is_symmetric() for M in B )
1836 Q
= QuaternionAlgebra(QQ
,-1,-1)
1839 # This is like the symmetric case, but we need to be careful:
1841 # * We want conjugate-symmetry, not just symmetry.
1842 # * The diagonal will (as a result) be real.
1846 for j
in range(i
+1):
1847 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1849 Sij
= cls
.real_embed(Eij
)
1852 # The second, third, and fourth ones have a minus
1853 # because they're conjugated.
1854 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1856 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1858 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1860 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1863 # Since we embedded these, we can drop back to the "field" that we
1864 # started with instead of the quaternion algebra "Q".
1865 return ( s
.change_ring(field
) for s
in S
)
1868 def __init__(self
, n
, field
=AA
, **kwargs
):
1869 basis
= self
._denormalized
_basis
(n
,field
)
1870 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1871 self
.rank
.set_cache(n
)
1874 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1876 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1877 with the half-trace inner product and jordan product ``x*y =
1878 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
1879 symmetric positive-definite "bilinear form" matrix. It has
1880 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
1881 when ``B`` is the identity matrix of order ``n-1``.
1885 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1886 ....: JordanSpinEJA)
1890 When no bilinear form is specified, the identity matrix is used,
1891 and the resulting algebra is the Jordan spin algebra::
1893 sage: J0 = BilinearFormEJA(3)
1894 sage: J1 = JordanSpinEJA(3)
1895 sage: J0.multiplication_table() == J0.multiplication_table()
1900 We can create a zero-dimensional algebra::
1902 sage: J = BilinearFormEJA(0)
1906 We can check the multiplication condition given in the Jordan, von
1907 Neumann, and Wigner paper (and also discussed on my "On the
1908 symmetry..." paper). Note that this relies heavily on the standard
1909 choice of basis, as does anything utilizing the bilinear form matrix::
1911 sage: set_random_seed()
1912 sage: n = ZZ.random_element(5)
1913 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
1914 sage: B = M.transpose()*M
1915 sage: J = BilinearFormEJA(n, B=B)
1916 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
1917 sage: V = J.vector_space()
1918 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
1919 ....: for ei in eis ]
1920 sage: actual = [ sis[i]*sis[j]
1921 ....: for i in range(n-1)
1922 ....: for j in range(n-1) ]
1923 sage: expected = [ J.one() if i == j else J.zero()
1924 ....: for i in range(n-1)
1925 ....: for j in range(n-1) ]
1926 sage: actual == expected
1929 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
1931 self
._B
= matrix
.identity(field
, max(0,n
-1))
1935 V
= VectorSpace(field
, n
)
1936 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1945 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
1946 zbar
= y0
*xbar
+ x0
*ybar
1947 z
= V([z0
] + zbar
.list())
1948 mult_table
[i
][j
] = z
1950 # The rank of this algebra is two, unless we're in a
1951 # one-dimensional ambient space (because the rank is bounded
1952 # by the ambient dimension).
1953 fdeja
= super(BilinearFormEJA
, self
)
1954 fdeja
.__init
__(field
, mult_table
, **kwargs
)
1955 self
.rank
.set_cache(min(n
,2))
1957 def inner_product(self
, x
, y
):
1959 Half of the trace inner product.
1961 This is defined so that the special case of the Jordan spin
1962 algebra gets the usual inner product.
1966 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1970 Ensure that this is one-half of the trace inner-product when
1971 the algebra isn't just the reals (when ``n`` isn't one). This
1972 is in Faraut and Koranyi, and also my "On the symmetry..."
1975 sage: set_random_seed()
1976 sage: n = ZZ.random_element(2,5)
1977 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
1978 sage: B = M.transpose()*M
1979 sage: J = BilinearFormEJA(n, B=B)
1980 sage: x = J.random_element()
1981 sage: y = J.random_element()
1982 sage: x.inner_product(y) == (x*y).trace()/2
1986 xvec
= x
.to_vector()
1988 yvec
= y
.to_vector()
1990 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
1993 class JordanSpinEJA(BilinearFormEJA
):
1995 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1996 with the usual inner product and jordan product ``x*y =
1997 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2002 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2006 This multiplication table can be verified by hand::
2008 sage: J = JordanSpinEJA(4)
2009 sage: e0,e1,e2,e3 = J.gens()
2025 We can change the generator prefix::
2027 sage: JordanSpinEJA(2, prefix='B').gens()
2032 Ensure that we have the usual inner product on `R^n`::
2034 sage: set_random_seed()
2035 sage: J = JordanSpinEJA.random_instance()
2036 sage: x,y = J.random_elements(2)
2037 sage: X = x.natural_representation()
2038 sage: Y = y.natural_representation()
2039 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2043 def __init__(self
, n
, field
=AA
, **kwargs
):
2044 # This is a special case of the BilinearFormEJA with the identity
2045 # matrix as its bilinear form.
2046 return super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2049 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2051 The trivial Euclidean Jordan algebra consisting of only a zero element.
2055 sage: from mjo.eja.eja_algebra import TrivialEJA
2059 sage: J = TrivialEJA()
2066 sage: 7*J.one()*12*J.one()
2068 sage: J.one().inner_product(J.one())
2070 sage: J.one().norm()
2072 sage: J.one().subalgebra_generated_by()
2073 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2078 def __init__(self
, field
=AA
, **kwargs
):
2080 fdeja
= super(TrivialEJA
, self
)
2081 # The rank is zero using my definition, namely the dimension of the
2082 # largest subalgebra generated by any element.
2083 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2084 self
.rank
.set_cache(0)