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.
10 sage: from mjo.eja.eja_algebra import random_eja
15 Euclidean Jordan algebra of dimension...
19 from itertools
import repeat
21 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
22 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
23 from sage
.combinat
.free_module
import CombinatorialFreeModule
24 from sage
.matrix
.constructor
import matrix
25 from sage
.matrix
.matrix_space
import MatrixSpace
26 from sage
.misc
.cachefunc
import cached_method
27 from sage
.misc
.lazy_import
import lazy_import
28 from sage
.misc
.table
import table
29 from sage
.modules
.free_module
import FreeModule
, VectorSpace
30 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
33 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
34 lazy_import('mjo.eja.eja_subalgebra',
35 'FiniteDimensionalEuclideanJordanSubalgebra')
36 from mjo
.eja
.eja_operator
import FiniteDimensionalEuclideanJordanAlgebraOperator
37 from mjo
.eja
.eja_utils
import _mat2vec
39 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
41 def _coerce_map_from_base_ring(self
):
43 Disable the map from the base ring into the algebra.
45 Performing a nonsense conversion like this automatically
46 is counterpedagogical. The fallback is to try the usual
47 element constructor, which should also fail.
51 sage: from mjo.eja.eja_algebra import random_eja
55 sage: set_random_seed()
56 sage: J = random_eja()
58 Traceback (most recent call last):
60 ValueError: not a naturally-represented algebra element
76 sage: from mjo.eja.eja_algebra import (
77 ....: FiniteDimensionalEuclideanJordanAlgebra,
83 By definition, Jordan multiplication commutes::
85 sage: set_random_seed()
86 sage: J = random_eja()
87 sage: x,y = J.random_elements(2)
93 The ``field`` we're given must be real with ``check_field=True``::
95 sage: JordanSpinEJA(2,QQbar)
96 Traceback (most recent call last):
98 ValueError: scalar field is not real
100 The multiplication table must be square with ``check_axioms=True``::
102 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
103 Traceback (most recent call last):
105 ValueError: multiplication table is not square
109 if not field
.is_subring(RR
):
110 # Note: this does return true for the real algebraic
111 # field, the rationals, and any quadratic field where
112 # we've specified a real embedding.
113 raise ValueError("scalar field is not real")
115 # The multiplication table had better be square
118 if not all( len(l
) == n
for l
in mult_table
):
119 raise ValueError("multiplication table is not square")
121 self
._natural
_basis
= natural_basis
124 category
= MagmaticAlgebras(field
).FiniteDimensional()
125 category
= category
.WithBasis().Unital()
127 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
132 self
.print_options(bracket
='')
134 # The multiplication table we're given is necessarily in terms
135 # of vectors, because we don't have an algebra yet for
136 # anything to be an element of. However, it's faster in the
137 # long run to have the multiplication table be in terms of
138 # algebra elements. We do this after calling the superclass
139 # constructor so that from_vector() knows what to do.
140 self
._multiplication
_table
= [
141 list(map(lambda x
: self
.from_vector(x
), ls
))
146 if not self
._is
_commutative
():
147 raise ValueError("algebra is not commutative")
148 if not self
._is
_jordanian
():
149 raise ValueError("Jordan identity does not hold")
150 if not self
._inner
_product
_is
_associative
():
151 raise ValueError("inner product is not associative")
153 def _element_constructor_(self
, elt
):
155 Construct an element of this algebra from its natural
158 This gets called only after the parent element _call_ method
159 fails to find a coercion for the argument.
163 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
165 ....: RealSymmetricEJA)
169 The identity in `S^n` is converted to the identity in the EJA::
171 sage: J = RealSymmetricEJA(3)
172 sage: I = matrix.identity(QQ,3)
173 sage: J(I) == J.one()
176 This skew-symmetric matrix can't be represented in the EJA::
178 sage: J = RealSymmetricEJA(3)
179 sage: A = matrix(QQ,3, lambda i,j: i-j)
181 Traceback (most recent call last):
183 ArithmeticError: vector is not in free module
187 Ensure that we can convert any element of the two non-matrix
188 simple algebras (whose natural representations are their usual
189 vector representations) back and forth faithfully::
191 sage: set_random_seed()
192 sage: J = HadamardEJA.random_instance()
193 sage: x = J.random_element()
194 sage: J(x.to_vector().column()) == x
196 sage: J = JordanSpinEJA.random_instance()
197 sage: x = J.random_element()
198 sage: J(x.to_vector().column()) == x
202 msg
= "not a naturally-represented algebra element"
204 # The superclass implementation of random_element()
205 # needs to be able to coerce "0" into the algebra.
207 elif elt
in self
.base_ring():
208 # Ensure that no base ring -> algebra coercion is performed
209 # by this method. There's some stupidity in sage that would
210 # otherwise propagate to this method; for example, sage thinks
211 # that the integer 3 belongs to the space of 2-by-2 matrices.
212 raise ValueError(msg
)
214 natural_basis
= self
.natural_basis()
215 basis_space
= natural_basis
[0].matrix_space()
216 if elt
not in basis_space
:
217 raise ValueError(msg
)
219 # Thanks for nothing! Matrix spaces aren't vector spaces in
220 # Sage, so we have to figure out its natural-basis coordinates
221 # ourselves. We use the basis space's ring instead of the
222 # element's ring because the basis space might be an algebraic
223 # closure whereas the base ring of the 3-by-3 identity matrix
224 # could be QQ instead of QQbar.
225 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
226 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
227 coords
= W
.coordinate_vector(_mat2vec(elt
))
228 return self
.from_vector(coords
)
232 Return a string representation of ``self``.
236 sage: from mjo.eja.eja_algebra import JordanSpinEJA
240 Ensure that it says what we think it says::
242 sage: JordanSpinEJA(2, field=AA)
243 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
244 sage: JordanSpinEJA(3, field=RDF)
245 Euclidean Jordan algebra of dimension 3 over Real Double Field
248 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
249 return fmt
.format(self
.dimension(), self
.base_ring())
251 def product_on_basis(self
, i
, j
):
252 return self
._multiplication
_table
[i
][j
]
254 def _is_commutative(self
):
256 Whether or not this algebra's multiplication table is commutative.
258 This method should of course always return ``True``, unless
259 this algebra was constructed with ``check_axioms=False`` and
260 passed an invalid multiplication table.
262 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
263 for i
in range(self
.dimension())
264 for j
in range(self
.dimension()) )
266 def _is_jordanian(self
):
268 Whether or not this algebra's multiplication table respects the
269 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
271 We only check one arrangement of `x` and `y`, so for a
272 ``True`` result to be truly true, you should also check
273 :meth:`_is_commutative`. This method should of course always
274 return ``True``, unless this algebra was constructed with
275 ``check_axioms=False`` and passed an invalid multiplication table.
277 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
279 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
280 for i
in range(self
.dimension())
281 for j
in range(self
.dimension()) )
283 def _inner_product_is_associative(self
):
285 Return whether or not this algebra's inner product `B` is
286 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
288 This method should of course always return ``True``, unless
289 this algebra was constructed with ``check_axioms=False`` and
290 passed an invalid multiplication table.
293 # Used to check whether or not something is zero in an inexact
294 # ring. This number is sufficient to allow the construction of
295 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
298 for i
in range(self
.dimension()):
299 for j
in range(self
.dimension()):
300 for k
in range(self
.dimension()):
304 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
306 if self
.base_ring().is_exact():
310 if diff
.abs() > epsilon
:
316 def characteristic_polynomial_of(self
):
318 Return the algebra's "characteristic polynomial of" function,
319 which is itself a multivariate polynomial that, when evaluated
320 at the coordinates of some algebra element, returns that
321 element's characteristic polynomial.
323 The resulting polynomial has `n+1` variables, where `n` is the
324 dimension of this algebra. The first `n` variables correspond to
325 the coordinates of an algebra element: when evaluated at the
326 coordinates of an algebra element with respect to a certain
327 basis, the result is a univariate polynomial (in the one
328 remaining variable ``t``), namely the characteristic polynomial
333 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
337 The characteristic polynomial in the spin algebra is given in
338 Alizadeh, Example 11.11::
340 sage: J = JordanSpinEJA(3)
341 sage: p = J.characteristic_polynomial_of(); p
342 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
343 sage: xvec = J.one().to_vector()
347 By definition, the characteristic polynomial is a monic
348 degree-zero polynomial in a rank-zero algebra. Note that
349 Cayley-Hamilton is indeed satisfied since the polynomial
350 ``1`` evaluates to the identity element of the algebra on
353 sage: J = TrivialEJA()
354 sage: J.characteristic_polynomial_of()
361 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
362 a
= self
._charpoly
_coefficients
()
364 # We go to a bit of trouble here to reorder the
365 # indeterminates, so that it's easier to evaluate the
366 # characteristic polynomial at x's coordinates and get back
367 # something in terms of t, which is what we want.
368 S
= PolynomialRing(self
.base_ring(),'t')
372 S
= PolynomialRing(S
, R
.variable_names())
375 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
378 def inner_product(self
, x
, y
):
380 The inner product associated with this Euclidean Jordan algebra.
382 Defaults to the trace inner product, but can be overridden by
383 subclasses if they are sure that the necessary properties are
388 sage: from mjo.eja.eja_algebra import random_eja
392 Our inner product is "associative," which means the following for
393 a symmetric bilinear form::
395 sage: set_random_seed()
396 sage: J = random_eja()
397 sage: x,y,z = J.random_elements(3)
398 sage: (x*y).inner_product(z) == y.inner_product(x*z)
402 X
= x
.natural_representation()
403 Y
= y
.natural_representation()
404 return self
.natural_inner_product(X
,Y
)
407 def is_trivial(self
):
409 Return whether or not this algebra is trivial.
411 A trivial algebra contains only the zero element.
415 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
420 sage: J = ComplexHermitianEJA(3)
426 sage: J = TrivialEJA()
431 return self
.dimension() == 0
434 def multiplication_table(self
):
436 Return a visual representation of this algebra's multiplication
437 table (on basis elements).
441 sage: from mjo.eja.eja_algebra import JordanSpinEJA
445 sage: J = JordanSpinEJA(4)
446 sage: J.multiplication_table()
447 +----++----+----+----+----+
448 | * || e0 | e1 | e2 | e3 |
449 +====++====+====+====+====+
450 | e0 || e0 | e1 | e2 | e3 |
451 +----++----+----+----+----+
452 | e1 || e1 | e0 | 0 | 0 |
453 +----++----+----+----+----+
454 | e2 || e2 | 0 | e0 | 0 |
455 +----++----+----+----+----+
456 | e3 || e3 | 0 | 0 | e0 |
457 +----++----+----+----+----+
460 M
= list(self
._multiplication
_table
) # copy
461 for i
in range(len(M
)):
462 # M had better be "square"
463 M
[i
] = [self
.monomial(i
)] + M
[i
]
464 M
= [["*"] + list(self
.gens())] + M
465 return table(M
, header_row
=True, header_column
=True, frame
=True)
468 def natural_basis(self
):
470 Return a more-natural representation of this algebra's basis.
472 Every finite-dimensional Euclidean Jordan Algebra is a direct
473 sum of five simple algebras, four of which comprise Hermitian
474 matrices. This method returns the original "natural" basis
475 for our underlying vector space. (Typically, the natural basis
476 is used to construct the multiplication table in the first place.)
478 Note that this will always return a matrix. The standard basis
479 in `R^n` will be returned as `n`-by-`1` column matrices.
483 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
484 ....: RealSymmetricEJA)
488 sage: J = RealSymmetricEJA(2)
490 Finite family {0: e0, 1: e1, 2: e2}
491 sage: J.natural_basis()
493 [1 0] [ 0 0.7071067811865475?] [0 0]
494 [0 0], [0.7071067811865475? 0], [0 1]
499 sage: J = JordanSpinEJA(2)
501 Finite family {0: e0, 1: e1}
502 sage: J.natural_basis()
509 if self
._natural
_basis
is None:
510 M
= self
.natural_basis_space()
511 return tuple( M(b
.to_vector()) for b
in self
.basis() )
513 return self
._natural
_basis
516 def natural_basis_space(self
):
518 Return the matrix space in which this algebra's natural basis
521 Generally this will be an `n`-by-`1` column-vector space,
522 except when the algebra is trivial. There it's `n`-by-`n`
523 (where `n` is zero), to ensure that two elements of the
524 natural basis space (empty matrices) can be multiplied.
526 if self
.is_trivial():
527 return MatrixSpace(self
.base_ring(), 0)
528 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
529 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
531 return self
._natural
_basis
[0].matrix_space()
535 def natural_inner_product(X
,Y
):
537 Compute the inner product of two naturally-represented elements.
539 For example in the real symmetric matrix EJA, this will compute
540 the trace inner-product of two n-by-n symmetric matrices. The
541 default should work for the real cartesian product EJA, the
542 Jordan spin EJA, and the real symmetric matrices. The others
543 will have to be overridden.
545 return (X
.conjugate_transpose()*Y
).trace()
551 Return the unit element of this algebra.
555 sage: from mjo.eja.eja_algebra import (HadamardEJA,
560 sage: J = HadamardEJA(5)
562 e0 + e1 + e2 + e3 + e4
566 The identity element acts like the identity::
568 sage: set_random_seed()
569 sage: J = random_eja()
570 sage: x = J.random_element()
571 sage: J.one()*x == x and x*J.one() == x
574 The matrix of the unit element's operator is the identity::
576 sage: set_random_seed()
577 sage: J = random_eja()
578 sage: actual = J.one().operator().matrix()
579 sage: expected = matrix.identity(J.base_ring(), J.dimension())
580 sage: actual == expected
583 Ensure that the cached unit element (often precomputed by
584 hand) agrees with the computed one::
586 sage: set_random_seed()
587 sage: J = random_eja()
588 sage: cached = J.one()
589 sage: J.one.clear_cache()
590 sage: J.one() == cached
594 # We can brute-force compute the matrices of the operators
595 # that correspond to the basis elements of this algebra.
596 # If some linear combination of those basis elements is the
597 # algebra identity, then the same linear combination of
598 # their matrices has to be the identity matrix.
600 # Of course, matrices aren't vectors in sage, so we have to
601 # appeal to the "long vectors" isometry.
602 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
604 # Now we use basic linear algebra to find the coefficients,
605 # of the matrices-as-vectors-linear-combination, which should
606 # work for the original algebra basis too.
607 A
= matrix(self
.base_ring(), oper_vecs
)
609 # We used the isometry on the left-hand side already, but we
610 # still need to do it for the right-hand side. Recall that we
611 # wanted something that summed to the identity matrix.
612 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
614 # Now if there's an identity element in the algebra, this
615 # should work. We solve on the left to avoid having to
616 # transpose the matrix "A".
617 return self
.from_vector(A
.solve_left(b
))
620 def peirce_decomposition(self
, c
):
622 The Peirce decomposition of this algebra relative to the
625 In the future, this can be extended to a complete system of
626 orthogonal idempotents.
630 - ``c`` -- an idempotent of this algebra.
634 A triple (J0, J5, J1) containing two subalgebras and one subspace
637 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
638 corresponding to the eigenvalue zero.
640 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
641 corresponding to the eigenvalue one-half.
643 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
644 corresponding to the eigenvalue one.
646 These are the only possible eigenspaces for that operator, and this
647 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
648 orthogonal, and are subalgebras of this algebra with the appropriate
653 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
657 The canonical example comes from the symmetric matrices, which
658 decompose into diagonal and off-diagonal parts::
660 sage: J = RealSymmetricEJA(3)
661 sage: C = matrix(QQ, [ [1,0,0],
665 sage: J0,J5,J1 = J.peirce_decomposition(c)
667 Euclidean Jordan algebra of dimension 1...
669 Vector space of degree 6 and dimension 2...
671 Euclidean Jordan algebra of dimension 3...
672 sage: J0.one().natural_representation()
676 sage: orig_df = AA.options.display_format
677 sage: AA.options.display_format = 'radical'
678 sage: J.from_vector(J5.basis()[0]).natural_representation()
682 sage: J.from_vector(J5.basis()[1]).natural_representation()
686 sage: AA.options.display_format = orig_df
687 sage: J1.one().natural_representation()
694 Every algebra decomposes trivially with respect to its identity
697 sage: set_random_seed()
698 sage: J = random_eja()
699 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
700 sage: J0.dimension() == 0 and J5.dimension() == 0
702 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
705 The decomposition is into eigenspaces, and its components are
706 therefore necessarily orthogonal. Moreover, the identity
707 elements in the two subalgebras are the projections onto their
708 respective subspaces of the superalgebra's identity element::
710 sage: set_random_seed()
711 sage: J = random_eja()
712 sage: x = J.random_element()
713 sage: if not J.is_trivial():
714 ....: while x.is_nilpotent():
715 ....: x = J.random_element()
716 sage: c = x.subalgebra_idempotent()
717 sage: J0,J5,J1 = J.peirce_decomposition(c)
719 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
720 ....: w = w.superalgebra_element()
721 ....: y = J.from_vector(y)
722 ....: z = z.superalgebra_element()
723 ....: ipsum += w.inner_product(y).abs()
724 ....: ipsum += w.inner_product(z).abs()
725 ....: ipsum += y.inner_product(z).abs()
728 sage: J1(c) == J1.one()
730 sage: J0(J.one() - c) == J0.one()
734 if not c
.is_idempotent():
735 raise ValueError("element is not idempotent: %s" % c
)
737 # Default these to what they should be if they turn out to be
738 # trivial, because eigenspaces_left() won't return eigenvalues
739 # corresponding to trivial spaces (e.g. it returns only the
740 # eigenspace corresponding to lambda=1 if you take the
741 # decomposition relative to the identity element).
742 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
743 J0
= trivial
# eigenvalue zero
744 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
745 J1
= trivial
# eigenvalue one
747 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
748 if eigval
== ~
(self
.base_ring()(2)):
751 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
752 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
760 raise ValueError("unexpected eigenvalue: %s" % eigval
)
765 def random_element(self
, thorough
=False):
767 Return a random element of this algebra.
769 Our algebra superclass method only returns a linear
770 combination of at most two basis elements. We instead
771 want the vector space "random element" method that
772 returns a more diverse selection.
776 - ``thorough`` -- (boolean; default False) whether or not we
777 should generate irrational coefficients for the random
778 element when our base ring is irrational; this slows the
779 algebra operations to a crawl, but any truly random method
783 # For a general base ring... maybe we can trust this to do the
784 # right thing? Unlikely, but.
785 V
= self
.vector_space()
786 v
= V
.random_element()
788 if self
.base_ring() is AA
:
789 # The "random element" method of the algebraic reals is
790 # stupid at the moment, and only returns integers between
791 # -2 and 2, inclusive:
793 # https://trac.sagemath.org/ticket/30875
795 # Instead, we implement our own "random vector" method,
796 # and then coerce that into the algebra. We use the vector
797 # space degree here instead of the dimension because a
798 # subalgebra could (for example) be spanned by only two
799 # vectors, each with five coordinates. We need to
800 # generate all five coordinates.
802 v
*= QQbar
.random_element().real()
804 v
*= QQ
.random_element()
806 return self
.from_vector(V
.coordinate_vector(v
))
808 def random_elements(self
, count
, thorough
=False):
810 Return ``count`` random elements as a tuple.
814 - ``thorough`` -- (boolean; default False) whether or not we
815 should generate irrational coefficients for the random
816 elements when our base ring is irrational; this slows the
817 algebra operations to a crawl, but any truly random method
822 sage: from mjo.eja.eja_algebra import JordanSpinEJA
826 sage: J = JordanSpinEJA(3)
827 sage: x,y,z = J.random_elements(3)
828 sage: all( [ x in J, y in J, z in J ])
830 sage: len( J.random_elements(10) ) == 10
834 return tuple( self
.random_element(thorough
)
835 for idx
in range(count
) )
839 def _charpoly_coefficients(self
):
841 The `r` polynomial coefficients of the "characteristic polynomial
845 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
846 R
= PolynomialRing(self
.base_ring(), var_names
)
848 F
= R
.fraction_field()
851 # From a result in my book, these are the entries of the
852 # basis representation of L_x.
853 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
856 L_x
= matrix(F
, n
, n
, L_x_i_j
)
859 if self
.rank
.is_in_cache():
861 # There's no need to pad the system with redundant
862 # columns if we *know* they'll be redundant.
865 # Compute an extra power in case the rank is equal to
866 # the dimension (otherwise, we would stop at x^(r-1)).
867 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
868 for k
in range(n
+1) ]
869 A
= matrix
.column(F
, x_powers
[:n
])
870 AE
= A
.extended_echelon_form()
877 # The theory says that only the first "r" coefficients are
878 # nonzero, and they actually live in the original polynomial
879 # ring and not the fraction field. We negate them because
880 # in the actual characteristic polynomial, they get moved
881 # to the other side where x^r lives.
882 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
887 Return the rank of this EJA.
889 This is a cached method because we know the rank a priori for
890 all of the algebras we can construct. Thus we can avoid the
891 expensive ``_charpoly_coefficients()`` call unless we truly
892 need to compute the whole characteristic polynomial.
896 sage: from mjo.eja.eja_algebra import (HadamardEJA,
898 ....: RealSymmetricEJA,
899 ....: ComplexHermitianEJA,
900 ....: QuaternionHermitianEJA,
905 The rank of the Jordan spin algebra is always two::
907 sage: JordanSpinEJA(2).rank()
909 sage: JordanSpinEJA(3).rank()
911 sage: JordanSpinEJA(4).rank()
914 The rank of the `n`-by-`n` Hermitian real, complex, or
915 quaternion matrices is `n`::
917 sage: RealSymmetricEJA(4).rank()
919 sage: ComplexHermitianEJA(3).rank()
921 sage: QuaternionHermitianEJA(2).rank()
926 Ensure that every EJA that we know how to construct has a
927 positive integer rank, unless the algebra is trivial in
928 which case its rank will be zero::
930 sage: set_random_seed()
931 sage: J = random_eja()
935 sage: r > 0 or (r == 0 and J.is_trivial())
938 Ensure that computing the rank actually works, since the ranks
939 of all simple algebras are known and will be cached by default::
941 sage: J = HadamardEJA(4)
942 sage: J.rank.clear_cache()
948 sage: J = JordanSpinEJA(4)
949 sage: J.rank.clear_cache()
955 sage: J = RealSymmetricEJA(3)
956 sage: J.rank.clear_cache()
962 sage: J = ComplexHermitianEJA(2)
963 sage: J.rank.clear_cache()
969 sage: J = QuaternionHermitianEJA(2)
970 sage: J.rank.clear_cache()
974 return len(self
._charpoly
_coefficients
())
977 def vector_space(self
):
979 Return the vector space that underlies this algebra.
983 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
987 sage: J = RealSymmetricEJA(2)
988 sage: J.vector_space()
989 Vector space of dimension 3 over...
992 return self
.zero().to_vector().parent().ambient_vector_space()
995 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
998 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1000 Algebras whose basis consists of vectors with rational
1001 entries. Equivalently, algebras whose multiplication tables
1002 contain only rational coefficients.
1004 When an EJA has a basis that can be made rational, we can speed up
1005 the computation of its characteristic polynomial by doing it over
1006 ``QQ``. All of the named EJA constructors that we provide fall
1010 def _charpoly_coefficients(self
):
1012 Override the parent method with something that tries to compute
1013 over a faster (non-extension) field.
1017 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1021 The base ring of the resulting polynomial coefficients is what
1022 it should be, and not the rationals (unless the algebra was
1023 already over the rationals)::
1025 sage: J = JordanSpinEJA(3)
1026 sage: J._charpoly_coefficients()
1027 (X1^2 - X2^2 - X3^2, -2*X1)
1028 sage: a0 = J._charpoly_coefficients()[0]
1030 Algebraic Real Field
1031 sage: a0.base_ring()
1032 Algebraic Real Field
1035 if self
.base_ring() is QQ
:
1036 # There's no need to construct *another* algebra over the
1037 # rationals if this one is already over the rationals.
1038 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1039 return superclass
._charpoly
_coefficients
()
1042 map(lambda x
: x
.to_vector(), ls
)
1043 for ls
in self
._multiplication
_table
1046 # Do the computation over the rationals. The answer will be
1047 # the same, because our basis coordinates are (essentially)
1049 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1053 a
= J
._charpoly
_coefficients
()
1054 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1057 class ConcreteEuclideanJordanAlgebra
:
1059 A class for the Euclidean Jordan algebras that we know by name.
1061 These are the Jordan algebras whose basis, multiplication table,
1062 rank, and so on are known a priori. More to the point, they are
1063 the Euclidean Jordan algebras for which we are able to conjure up
1064 a "random instance."
1068 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1072 Our natural basis is normalized with respect to the natural inner
1073 product unless we specify otherwise::
1075 sage: set_random_seed()
1076 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1077 sage: all( b.norm() == 1 for b in J.gens() )
1080 Since our natural basis is normalized with respect to the natural
1081 inner product, and since we know that this algebra is an EJA, any
1082 left-multiplication operator's matrix will be symmetric because
1083 natural->EJA basis representation is an isometry and within the EJA
1084 the operator is self-adjoint by the Jordan axiom::
1086 sage: set_random_seed()
1087 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1088 sage: x = J.random_element()
1089 sage: x.operator().matrix().is_symmetric()
1095 def _max_random_instance_size():
1097 Return an integer "size" that is an upper bound on the size of
1098 this algebra when it is used in a random test
1099 case. Unfortunately, the term "size" is ambiguous -- when
1100 dealing with `R^n` under either the Hadamard or Jordan spin
1101 product, the "size" refers to the dimension `n`. When dealing
1102 with a matrix algebra (real symmetric or complex/quaternion
1103 Hermitian), it refers to the size of the matrix, which is far
1104 less than the dimension of the underlying vector space.
1106 This method must be implemented in each subclass.
1108 raise NotImplementedError
1111 def random_instance(cls
, field
=AA
, **kwargs
):
1113 Return a random instance of this type of algebra.
1115 This method should be implemented in each subclass.
1117 from sage
.misc
.prandom
import choice
1118 eja_class
= choice(cls
.__subclasses
__())
1119 return eja_class
.random_instance(field
)
1122 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1124 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1126 Compared to the superclass constructor, we take a basis instead of
1127 a multiplication table because the latter can be computed in terms
1128 of the former when the product is known (like it is here).
1130 # Used in this class's fast _charpoly_coefficients() override.
1131 self
._basis
_normalizers
= None
1133 # We're going to loop through this a few times, so now's a good
1134 # time to ensure that it isn't a generator expression.
1135 basis
= tuple(basis
)
1137 algebra_dim
= len(basis
)
1138 if algebra_dim
> 1 and normalize_basis
:
1139 # We'll need sqrt(2) to normalize the basis, and this
1140 # winds up in the multiplication table, so the whole
1141 # algebra needs to be over the field extension.
1142 R
= PolynomialRing(field
, 'z')
1145 if p
.is_irreducible():
1146 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1147 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1148 self
._basis
_normalizers
= tuple(
1149 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1150 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1152 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1154 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1156 natural_basis
=basis
,
1159 if algebra_dim
== 0:
1160 self
.one
.set_cache(self
.zero())
1162 n
= basis
[0].nrows()
1163 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1164 # details of this algebra.
1165 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1169 def _charpoly_coefficients(self
):
1171 Override the parent method with something that tries to compute
1172 over a faster (non-extension) field.
1174 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1175 # We didn't normalize, or the basis we started with had
1176 # entries in a nice field already. Just compute the thing.
1177 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1179 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1180 self
._basis
_normalizers
) )
1182 # Do this over the rationals and convert back at the end.
1183 # Only works because we know the entries of the basis are
1184 # integers. The argument ``check_axioms=False`` is required
1185 # because the trace inner-product method for this
1186 # class is a stub and can't actually be checked.
1187 J
= MatrixEuclideanJordanAlgebra(QQ
,
1189 normalize_basis
=False,
1192 a
= J
._charpoly
_coefficients
()
1194 # Unfortunately, changing the basis does change the
1195 # coefficients of the characteristic polynomial, but since
1196 # these are really the coefficients of the "characteristic
1197 # polynomial of" function, everything is still nice and
1198 # unevaluated. It's therefore "obvious" how scaling the
1199 # basis affects the coordinate variables X1, X2, et
1200 # cetera. Scaling the first basis vector up by "n" adds a
1201 # factor of 1/n into every "X1" term, for example. So here
1202 # we simply undo the basis_normalizer scaling that we
1203 # performed earlier.
1205 # The a[0] access here is safe because trivial algebras
1206 # won't have any basis normalizers and therefore won't
1207 # make it to this "else" branch.
1208 XS
= a
[0].parent().gens()
1209 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1210 for i
in range(len(XS
)) }
1211 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1215 def multiplication_table_from_matrix_basis(basis
):
1217 At least three of the five simple Euclidean Jordan algebras have the
1218 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1219 multiplication on the right is matrix multiplication. Given a basis
1220 for the underlying matrix space, this function returns a
1221 multiplication table (obtained by looping through the basis
1222 elements) for an algebra of those matrices.
1224 # In S^2, for example, we nominally have four coordinates even
1225 # though the space is of dimension three only. The vector space V
1226 # is supposed to hold the entire long vector, and the subspace W
1227 # of V will be spanned by the vectors that arise from symmetric
1228 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1232 field
= basis
[0].base_ring()
1233 dimension
= basis
[0].nrows()
1235 V
= VectorSpace(field
, dimension
**2)
1236 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1238 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1241 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1242 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1250 Embed the matrix ``M`` into a space of real matrices.
1252 The matrix ``M`` can have entries in any field at the moment:
1253 the real numbers, complex numbers, or quaternions. And although
1254 they are not a field, we can probably support octonions at some
1255 point, too. This function returns a real matrix that "acts like"
1256 the original with respect to matrix multiplication; i.e.
1258 real_embed(M*N) = real_embed(M)*real_embed(N)
1261 raise NotImplementedError
1265 def real_unembed(M
):
1267 The inverse of :meth:`real_embed`.
1269 raise NotImplementedError
1273 def natural_inner_product(cls
,X
,Y
):
1274 Xu
= cls
.real_unembed(X
)
1275 Yu
= cls
.real_unembed(Y
)
1276 tr
= (Xu
*Yu
).trace()
1279 # Works in QQ, AA, RDF, et cetera.
1281 except AttributeError:
1282 # A quaternion doesn't have a real() method, but does
1283 # have coefficient_tuple() method that returns the
1284 # coefficients of 1, i, j, and k -- in that order.
1285 return tr
.coefficient_tuple()[0]
1288 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1292 The identity function, for embedding real matrices into real
1298 def real_unembed(M
):
1300 The identity function, for unembedding real matrices from real
1306 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1307 ConcreteEuclideanJordanAlgebra
):
1309 The rank-n simple EJA consisting of real symmetric n-by-n
1310 matrices, the usual symmetric Jordan product, and the trace inner
1311 product. It has dimension `(n^2 + n)/2` over the reals.
1315 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1319 sage: J = RealSymmetricEJA(2)
1320 sage: e0, e1, e2 = J.gens()
1328 In theory, our "field" can be any subfield of the reals::
1330 sage: RealSymmetricEJA(2, RDF)
1331 Euclidean Jordan algebra of dimension 3 over Real Double Field
1332 sage: RealSymmetricEJA(2, RR)
1333 Euclidean Jordan algebra of dimension 3 over Real Field with
1334 53 bits of precision
1338 The dimension of this algebra is `(n^2 + n) / 2`::
1340 sage: set_random_seed()
1341 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1342 sage: n = ZZ.random_element(1, n_max)
1343 sage: J = RealSymmetricEJA(n)
1344 sage: J.dimension() == (n^2 + n)/2
1347 The Jordan multiplication is what we think it is::
1349 sage: set_random_seed()
1350 sage: J = RealSymmetricEJA.random_instance()
1351 sage: x,y = J.random_elements(2)
1352 sage: actual = (x*y).natural_representation()
1353 sage: X = x.natural_representation()
1354 sage: Y = y.natural_representation()
1355 sage: expected = (X*Y + Y*X)/2
1356 sage: actual == expected
1358 sage: J(expected) == x*y
1361 We can change the generator prefix::
1363 sage: RealSymmetricEJA(3, prefix='q').gens()
1364 (q0, q1, q2, q3, q4, q5)
1366 We can construct the (trivial) algebra of rank zero::
1368 sage: RealSymmetricEJA(0)
1369 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1373 def _denormalized_basis(cls
, n
, field
):
1375 Return a basis for the space of real symmetric n-by-n matrices.
1379 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1383 sage: set_random_seed()
1384 sage: n = ZZ.random_element(1,5)
1385 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1386 sage: all( M.is_symmetric() for M in B)
1390 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1394 for j
in range(i
+1):
1395 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1399 Sij
= Eij
+ Eij
.transpose()
1405 def _max_random_instance_size():
1406 return 4 # Dimension 10
1409 def random_instance(cls
, field
=AA
, **kwargs
):
1411 Return a random instance of this type of algebra.
1413 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1414 return cls(n
, field
, **kwargs
)
1416 def __init__(self
, n
, field
=AA
, **kwargs
):
1417 basis
= self
._denormalized
_basis
(n
, field
)
1418 super(RealSymmetricEJA
, self
).__init
__(field
,
1422 self
.rank
.set_cache(n
)
1425 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1429 Embed the n-by-n complex matrix ``M`` into the space of real
1430 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1431 bi` to the block matrix ``[[a,b],[-b,a]]``.
1435 sage: from mjo.eja.eja_algebra import \
1436 ....: ComplexMatrixEuclideanJordanAlgebra
1440 sage: F = QuadraticField(-1, 'I')
1441 sage: x1 = F(4 - 2*i)
1442 sage: x2 = F(1 + 2*i)
1445 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1446 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1455 Embedding is a homomorphism (isomorphism, in fact)::
1457 sage: set_random_seed()
1458 sage: n = ZZ.random_element(3)
1459 sage: F = QuadraticField(-1, 'I')
1460 sage: X = random_matrix(F, n)
1461 sage: Y = random_matrix(F, n)
1462 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1463 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1464 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1471 raise ValueError("the matrix 'M' must be square")
1473 # We don't need any adjoined elements...
1474 field
= M
.base_ring().base_ring()
1478 a
= z
.list()[0] # real part, I guess
1479 b
= z
.list()[1] # imag part, I guess
1480 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1482 return matrix
.block(field
, n
, blocks
)
1486 def real_unembed(M
):
1488 The inverse of _embed_complex_matrix().
1492 sage: from mjo.eja.eja_algebra import \
1493 ....: ComplexMatrixEuclideanJordanAlgebra
1497 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1498 ....: [-2, 1, -4, 3],
1499 ....: [ 9, 10, 11, 12],
1500 ....: [-10, 9, -12, 11] ])
1501 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1503 [ 10*I + 9 12*I + 11]
1507 Unembedding is the inverse of embedding::
1509 sage: set_random_seed()
1510 sage: F = QuadraticField(-1, 'I')
1511 sage: M = random_matrix(F, 3)
1512 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1513 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1519 raise ValueError("the matrix 'M' must be square")
1520 if not n
.mod(2).is_zero():
1521 raise ValueError("the matrix 'M' must be a complex embedding")
1523 # If "M" was normalized, its base ring might have roots
1524 # adjoined and they can stick around after unembedding.
1525 field
= M
.base_ring()
1526 R
= PolynomialRing(field
, 'z')
1529 # Sage doesn't know how to embed AA into QQbar, i.e. how
1530 # to adjoin sqrt(-1) to AA.
1533 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1536 # Go top-left to bottom-right (reading order), converting every
1537 # 2-by-2 block we see to a single complex element.
1539 for k
in range(n
/2):
1540 for j
in range(n
/2):
1541 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1542 if submat
[0,0] != submat
[1,1]:
1543 raise ValueError('bad on-diagonal submatrix')
1544 if submat
[0,1] != -submat
[1,0]:
1545 raise ValueError('bad off-diagonal submatrix')
1546 z
= submat
[0,0] + submat
[0,1]*i
1549 return matrix(F
, n
/2, elements
)
1553 def natural_inner_product(cls
,X
,Y
):
1555 Compute a natural inner product in this algebra directly from
1560 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1564 This gives the same answer as the slow, default method implemented
1565 in :class:`MatrixEuclideanJordanAlgebra`::
1567 sage: set_random_seed()
1568 sage: J = ComplexHermitianEJA.random_instance()
1569 sage: x,y = J.random_elements(2)
1570 sage: Xe = x.natural_representation()
1571 sage: Ye = y.natural_representation()
1572 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1573 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1574 sage: expected = (X*Y).trace().real()
1575 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1576 sage: actual == expected
1580 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1583 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1584 ConcreteEuclideanJordanAlgebra
):
1586 The rank-n simple EJA consisting of complex Hermitian n-by-n
1587 matrices over the real numbers, the usual symmetric Jordan product,
1588 and the real-part-of-trace inner product. It has dimension `n^2` over
1593 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1597 In theory, our "field" can be any subfield of the reals::
1599 sage: ComplexHermitianEJA(2, RDF)
1600 Euclidean Jordan algebra of dimension 4 over Real Double Field
1601 sage: ComplexHermitianEJA(2, RR)
1602 Euclidean Jordan algebra of dimension 4 over Real Field with
1603 53 bits of precision
1607 The dimension of this algebra is `n^2`::
1609 sage: set_random_seed()
1610 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1611 sage: n = ZZ.random_element(1, n_max)
1612 sage: J = ComplexHermitianEJA(n)
1613 sage: J.dimension() == n^2
1616 The Jordan multiplication is what we think it is::
1618 sage: set_random_seed()
1619 sage: J = ComplexHermitianEJA.random_instance()
1620 sage: x,y = J.random_elements(2)
1621 sage: actual = (x*y).natural_representation()
1622 sage: X = x.natural_representation()
1623 sage: Y = y.natural_representation()
1624 sage: expected = (X*Y + Y*X)/2
1625 sage: actual == expected
1627 sage: J(expected) == x*y
1630 We can change the generator prefix::
1632 sage: ComplexHermitianEJA(2, prefix='z').gens()
1635 We can construct the (trivial) algebra of rank zero::
1637 sage: ComplexHermitianEJA(0)
1638 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1643 def _denormalized_basis(cls
, n
, field
):
1645 Returns a basis for the space of complex Hermitian n-by-n matrices.
1647 Why do we embed these? Basically, because all of numerical linear
1648 algebra assumes that you're working with vectors consisting of `n`
1649 entries from a field and scalars from the same field. There's no way
1650 to tell SageMath that (for example) the vectors contain complex
1651 numbers, while the scalar field is real.
1655 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1659 sage: set_random_seed()
1660 sage: n = ZZ.random_element(1,5)
1661 sage: field = QuadraticField(2, 'sqrt2')
1662 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1663 sage: all( M.is_symmetric() for M in B)
1667 R
= PolynomialRing(field
, 'z')
1669 F
= field
.extension(z
**2 + 1, 'I')
1672 # This is like the symmetric case, but we need to be careful:
1674 # * We want conjugate-symmetry, not just symmetry.
1675 # * The diagonal will (as a result) be real.
1679 for j
in range(i
+1):
1680 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1682 Sij
= cls
.real_embed(Eij
)
1685 # The second one has a minus because it's conjugated.
1686 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1688 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1691 # Since we embedded these, we can drop back to the "field" that we
1692 # started with instead of the complex extension "F".
1693 return ( s
.change_ring(field
) for s
in S
)
1696 def __init__(self
, n
, field
=AA
, **kwargs
):
1697 basis
= self
._denormalized
_basis
(n
,field
)
1698 super(ComplexHermitianEJA
,self
).__init
__(field
,
1702 self
.rank
.set_cache(n
)
1705 def _max_random_instance_size():
1706 return 3 # Dimension 9
1709 def random_instance(cls
, field
=AA
, **kwargs
):
1711 Return a random instance of this type of algebra.
1713 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1714 return cls(n
, field
, **kwargs
)
1716 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1720 Embed the n-by-n quaternion matrix ``M`` into the space of real
1721 matrices of size 4n-by-4n by first sending each quaternion entry `z
1722 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1723 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1728 sage: from mjo.eja.eja_algebra import \
1729 ....: QuaternionMatrixEuclideanJordanAlgebra
1733 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1734 sage: i,j,k = Q.gens()
1735 sage: x = 1 + 2*i + 3*j + 4*k
1736 sage: M = matrix(Q, 1, [[x]])
1737 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1743 Embedding is a homomorphism (isomorphism, in fact)::
1745 sage: set_random_seed()
1746 sage: n = ZZ.random_element(2)
1747 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1748 sage: X = random_matrix(Q, n)
1749 sage: Y = random_matrix(Q, n)
1750 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1751 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1752 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1757 quaternions
= M
.base_ring()
1760 raise ValueError("the matrix 'M' must be square")
1762 F
= QuadraticField(-1, 'I')
1767 t
= z
.coefficient_tuple()
1772 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1773 [-c
+ d
*i
, a
- b
*i
]])
1774 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1775 blocks
.append(realM
)
1777 # We should have real entries by now, so use the realest field
1778 # we've got for the return value.
1779 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1784 def real_unembed(M
):
1786 The inverse of _embed_quaternion_matrix().
1790 sage: from mjo.eja.eja_algebra import \
1791 ....: QuaternionMatrixEuclideanJordanAlgebra
1795 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1796 ....: [-2, 1, -4, 3],
1797 ....: [-3, 4, 1, -2],
1798 ....: [-4, -3, 2, 1]])
1799 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1800 [1 + 2*i + 3*j + 4*k]
1804 Unembedding is the inverse of embedding::
1806 sage: set_random_seed()
1807 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1808 sage: M = random_matrix(Q, 3)
1809 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1810 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1816 raise ValueError("the matrix 'M' must be square")
1817 if not n
.mod(4).is_zero():
1818 raise ValueError("the matrix 'M' must be a quaternion embedding")
1820 # Use the base ring of the matrix to ensure that its entries can be
1821 # multiplied by elements of the quaternion algebra.
1822 field
= M
.base_ring()
1823 Q
= QuaternionAlgebra(field
,-1,-1)
1826 # Go top-left to bottom-right (reading order), converting every
1827 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1830 for l
in range(n
/4):
1831 for m
in range(n
/4):
1832 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1833 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1834 if submat
[0,0] != submat
[1,1].conjugate():
1835 raise ValueError('bad on-diagonal submatrix')
1836 if submat
[0,1] != -submat
[1,0].conjugate():
1837 raise ValueError('bad off-diagonal submatrix')
1838 z
= submat
[0,0].real()
1839 z
+= submat
[0,0].imag()*i
1840 z
+= submat
[0,1].real()*j
1841 z
+= submat
[0,1].imag()*k
1844 return matrix(Q
, n
/4, elements
)
1848 def natural_inner_product(cls
,X
,Y
):
1850 Compute a natural inner product in this algebra directly from
1855 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1859 This gives the same answer as the slow, default method implemented
1860 in :class:`MatrixEuclideanJordanAlgebra`::
1862 sage: set_random_seed()
1863 sage: J = QuaternionHermitianEJA.random_instance()
1864 sage: x,y = J.random_elements(2)
1865 sage: Xe = x.natural_representation()
1866 sage: Ye = y.natural_representation()
1867 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1868 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1869 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1870 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1871 sage: actual == expected
1875 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1878 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1879 ConcreteEuclideanJordanAlgebra
):
1881 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1882 matrices, the usual symmetric Jordan product, and the
1883 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1888 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1892 In theory, our "field" can be any subfield of the reals::
1894 sage: QuaternionHermitianEJA(2, RDF)
1895 Euclidean Jordan algebra of dimension 6 over Real Double Field
1896 sage: QuaternionHermitianEJA(2, RR)
1897 Euclidean Jordan algebra of dimension 6 over Real Field with
1898 53 bits of precision
1902 The dimension of this algebra is `2*n^2 - n`::
1904 sage: set_random_seed()
1905 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
1906 sage: n = ZZ.random_element(1, n_max)
1907 sage: J = QuaternionHermitianEJA(n)
1908 sage: J.dimension() == 2*(n^2) - n
1911 The Jordan multiplication is what we think it is::
1913 sage: set_random_seed()
1914 sage: J = QuaternionHermitianEJA.random_instance()
1915 sage: x,y = J.random_elements(2)
1916 sage: actual = (x*y).natural_representation()
1917 sage: X = x.natural_representation()
1918 sage: Y = y.natural_representation()
1919 sage: expected = (X*Y + Y*X)/2
1920 sage: actual == expected
1922 sage: J(expected) == x*y
1925 We can change the generator prefix::
1927 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1928 (a0, a1, a2, a3, a4, a5)
1930 We can construct the (trivial) algebra of rank zero::
1932 sage: QuaternionHermitianEJA(0)
1933 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1937 def _denormalized_basis(cls
, n
, field
):
1939 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1941 Why do we embed these? Basically, because all of numerical
1942 linear algebra assumes that you're working with vectors consisting
1943 of `n` entries from a field and scalars from the same field. There's
1944 no way to tell SageMath that (for example) the vectors contain
1945 complex numbers, while the scalar field is real.
1949 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1953 sage: set_random_seed()
1954 sage: n = ZZ.random_element(1,5)
1955 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1956 sage: all( M.is_symmetric() for M in B )
1960 Q
= QuaternionAlgebra(QQ
,-1,-1)
1963 # This is like the symmetric case, but we need to be careful:
1965 # * We want conjugate-symmetry, not just symmetry.
1966 # * The diagonal will (as a result) be real.
1970 for j
in range(i
+1):
1971 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1973 Sij
= cls
.real_embed(Eij
)
1976 # The second, third, and fourth ones have a minus
1977 # because they're conjugated.
1978 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1980 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1982 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1984 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1987 # Since we embedded these, we can drop back to the "field" that we
1988 # started with instead of the quaternion algebra "Q".
1989 return ( s
.change_ring(field
) for s
in S
)
1992 def __init__(self
, n
, field
=AA
, **kwargs
):
1993 basis
= self
._denormalized
_basis
(n
,field
)
1994 super(QuaternionHermitianEJA
,self
).__init
__(field
,
1998 self
.rank
.set_cache(n
)
2001 def _max_random_instance_size():
2003 The maximum rank of a random QuaternionHermitianEJA.
2005 return 2 # Dimension 6
2008 def random_instance(cls
, field
=AA
, **kwargs
):
2010 Return a random instance of this type of algebra.
2012 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2013 return cls(n
, field
, **kwargs
)
2016 class HadamardEJA(RationalBasisEuclideanJordanAlgebra
,
2017 ConcreteEuclideanJordanAlgebra
):
2019 Return the Euclidean Jordan Algebra corresponding to the set
2020 `R^n` under the Hadamard product.
2022 Note: this is nothing more than the Cartesian product of ``n``
2023 copies of the spin algebra. Once Cartesian product algebras
2024 are implemented, this can go.
2028 sage: from mjo.eja.eja_algebra import HadamardEJA
2032 This multiplication table can be verified by hand::
2034 sage: J = HadamardEJA(3)
2035 sage: e0,e1,e2 = J.gens()
2051 We can change the generator prefix::
2053 sage: HadamardEJA(3, prefix='r').gens()
2057 def __init__(self
, n
, field
=AA
, **kwargs
):
2058 V
= VectorSpace(field
, n
)
2059 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
2062 super(HadamardEJA
, self
).__init
__(field
,
2066 self
.rank
.set_cache(n
)
2069 self
.one
.set_cache( self
.zero() )
2071 self
.one
.set_cache( sum(self
.gens()) )
2074 def _max_random_instance_size():
2076 The maximum dimension of a random HadamardEJA.
2081 def random_instance(cls
, field
=AA
, **kwargs
):
2083 Return a random instance of this type of algebra.
2085 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2086 return cls(n
, field
, **kwargs
)
2089 def inner_product(self
, x
, y
):
2091 Faster to reimplement than to use natural representations.
2095 sage: from mjo.eja.eja_algebra import HadamardEJA
2099 Ensure that this is the usual inner product for the algebras
2102 sage: set_random_seed()
2103 sage: J = HadamardEJA.random_instance()
2104 sage: x,y = J.random_elements(2)
2105 sage: X = x.natural_representation()
2106 sage: Y = y.natural_representation()
2107 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2111 return x
.to_vector().inner_product(y
.to_vector())
2114 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
,
2115 ConcreteEuclideanJordanAlgebra
):
2117 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2118 with the half-trace inner product and jordan product ``x*y =
2119 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2120 a symmetric positive-definite "bilinear form" matrix. Its
2121 dimension is the size of `B`, and it has rank two in dimensions
2122 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2123 the identity matrix of order ``n``.
2125 We insist that the one-by-one upper-left identity block of `B` be
2126 passed in as well so that we can be passed a matrix of size zero
2127 to construct a trivial algebra.
2131 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2132 ....: JordanSpinEJA)
2136 When no bilinear form is specified, the identity matrix is used,
2137 and the resulting algebra is the Jordan spin algebra::
2139 sage: B = matrix.identity(AA,3)
2140 sage: J0 = BilinearFormEJA(B)
2141 sage: J1 = JordanSpinEJA(3)
2142 sage: J0.multiplication_table() == J0.multiplication_table()
2145 An error is raised if the matrix `B` does not correspond to a
2146 positive-definite bilinear form::
2148 sage: B = matrix.random(QQ,2,3)
2149 sage: J = BilinearFormEJA(B)
2150 Traceback (most recent call last):
2152 ValueError: bilinear form is not positive-definite
2153 sage: B = matrix.zero(QQ,3)
2154 sage: J = BilinearFormEJA(B)
2155 Traceback (most recent call last):
2157 ValueError: bilinear form is not positive-definite
2161 We can create a zero-dimensional algebra::
2163 sage: B = matrix.identity(AA,0)
2164 sage: J = BilinearFormEJA(B)
2168 We can check the multiplication condition given in the Jordan, von
2169 Neumann, and Wigner paper (and also discussed on my "On the
2170 symmetry..." paper). Note that this relies heavily on the standard
2171 choice of basis, as does anything utilizing the bilinear form matrix::
2173 sage: set_random_seed()
2174 sage: n = ZZ.random_element(5)
2175 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2176 sage: B11 = matrix.identity(QQ,1)
2177 sage: B22 = M.transpose()*M
2178 sage: B = block_matrix(2,2,[ [B11,0 ],
2180 sage: J = BilinearFormEJA(B)
2181 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2182 sage: V = J.vector_space()
2183 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2184 ....: for ei in eis ]
2185 sage: actual = [ sis[i]*sis[j]
2186 ....: for i in range(n-1)
2187 ....: for j in range(n-1) ]
2188 sage: expected = [ J.one() if i == j else J.zero()
2189 ....: for i in range(n-1)
2190 ....: for j in range(n-1) ]
2191 sage: actual == expected
2194 def __init__(self
, B
, field
=AA
, **kwargs
):
2198 if not B
.is_positive_definite():
2199 raise ValueError("bilinear form is not positive-definite")
2201 V
= VectorSpace(field
, n
)
2202 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2211 z0
= (B
*x
).inner_product(y
)
2212 zbar
= y0
*xbar
+ x0
*ybar
2213 z
= V([z0
] + zbar
.list())
2214 mult_table
[i
][j
] = z
2216 # The rank of this algebra is two, unless we're in a
2217 # one-dimensional ambient space (because the rank is bounded
2218 # by the ambient dimension).
2219 super(BilinearFormEJA
, self
).__init
__(field
,
2223 self
.rank
.set_cache(min(n
,2))
2226 self
.one
.set_cache( self
.zero() )
2228 self
.one
.set_cache( self
.monomial(0) )
2231 def _max_random_instance_size():
2233 The maximum dimension of a random BilinearFormEJA.
2238 def random_instance(cls
, field
=AA
, **kwargs
):
2240 Return a random instance of this algebra.
2242 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2244 B
= matrix
.identity(field
, n
)
2245 return cls(B
, field
, **kwargs
)
2247 B11
= matrix
.identity(field
,1)
2248 M
= matrix
.random(field
, n
-1)
2249 I
= matrix
.identity(field
, n
-1)
2250 alpha
= field
.zero()
2251 while alpha
.is_zero():
2252 alpha
= field
.random_element().abs()
2253 B22
= M
.transpose()*M
+ alpha
*I
2255 from sage
.matrix
.special
import block_matrix
2256 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2259 return cls(B
, field
, **kwargs
)
2261 def inner_product(self
, x
, y
):
2263 Half of the trace inner product.
2265 This is defined so that the special case of the Jordan spin
2266 algebra gets the usual inner product.
2270 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2274 Ensure that this is one-half of the trace inner-product when
2275 the algebra isn't just the reals (when ``n`` isn't one). This
2276 is in Faraut and Koranyi, and also my "On the symmetry..."
2279 sage: set_random_seed()
2280 sage: J = BilinearFormEJA.random_instance()
2281 sage: n = J.dimension()
2282 sage: x = J.random_element()
2283 sage: y = J.random_element()
2284 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
2288 return (self
._B
*x
.to_vector()).inner_product(y
.to_vector())
2291 class JordanSpinEJA(BilinearFormEJA
):
2293 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2294 with the usual inner product and jordan product ``x*y =
2295 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2300 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2304 This multiplication table can be verified by hand::
2306 sage: J = JordanSpinEJA(4)
2307 sage: e0,e1,e2,e3 = J.gens()
2323 We can change the generator prefix::
2325 sage: JordanSpinEJA(2, prefix='B').gens()
2330 Ensure that we have the usual inner product on `R^n`::
2332 sage: set_random_seed()
2333 sage: J = JordanSpinEJA.random_instance()
2334 sage: x,y = J.random_elements(2)
2335 sage: X = x.natural_representation()
2336 sage: Y = y.natural_representation()
2337 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2341 def __init__(self
, n
, field
=AA
, **kwargs
):
2342 # This is a special case of the BilinearFormEJA with the identity
2343 # matrix as its bilinear form.
2344 B
= matrix
.identity(field
, n
)
2345 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2348 def _max_random_instance_size():
2350 The maximum dimension of a random JordanSpinEJA.
2355 def random_instance(cls
, field
=AA
, **kwargs
):
2357 Return a random instance of this type of algebra.
2359 Needed here to override the implementation for ``BilinearFormEJA``.
2361 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2362 return cls(n
, field
, **kwargs
)
2365 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
,
2366 ConcreteEuclideanJordanAlgebra
):
2368 The trivial Euclidean Jordan algebra consisting of only a zero element.
2372 sage: from mjo.eja.eja_algebra import TrivialEJA
2376 sage: J = TrivialEJA()
2383 sage: 7*J.one()*12*J.one()
2385 sage: J.one().inner_product(J.one())
2387 sage: J.one().norm()
2389 sage: J.one().subalgebra_generated_by()
2390 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2395 def __init__(self
, field
=AA
, **kwargs
):
2397 super(TrivialEJA
, self
).__init
__(field
,
2401 # The rank is zero using my definition, namely the dimension of the
2402 # largest subalgebra generated by any element.
2403 self
.rank
.set_cache(0)
2404 self
.one
.set_cache( self
.zero() )
2407 def random_instance(cls
, field
=AA
, **kwargs
):
2408 # We don't take a "size" argument so the superclass method is
2409 # inappropriate for us.
2410 return cls(field
, **kwargs
)
2412 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2414 The external (orthogonal) direct sum of two other Euclidean Jordan
2415 algebras. Essentially the Cartesian product of its two factors.
2416 Every Euclidean Jordan algebra decomposes into an orthogonal
2417 direct sum of simple Euclidean Jordan algebras, so no generality
2418 is lost by providing only this construction.
2422 sage: from mjo.eja.eja_algebra import (random_eja,
2424 ....: RealSymmetricEJA,
2429 sage: J1 = HadamardEJA(2)
2430 sage: J2 = RealSymmetricEJA(3)
2431 sage: J = DirectSumEJA(J1,J2)
2439 The external direct sum construction is only valid when the two factors
2440 have the same base ring; an error is raised otherwise::
2442 sage: set_random_seed()
2443 sage: J1 = random_eja(AA)
2444 sage: J2 = random_eja(QQ)
2445 sage: J = DirectSumEJA(J1,J2)
2446 Traceback (most recent call last):
2448 ValueError: algebras must share the same base field
2451 def __init__(self
, J1
, J2
, **kwargs
):
2452 if J1
.base_ring() != J2
.base_ring():
2453 raise ValueError("algebras must share the same base field")
2454 field
= J1
.base_ring()
2456 self
._factors
= (J1
, J2
)
2460 V
= VectorSpace(field
, n
)
2461 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2465 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2466 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2470 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2471 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2473 super(DirectSumEJA
, self
).__init
__(field
,
2477 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2482 Return the pair of this algebra's factors.
2486 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2487 ....: JordanSpinEJA,
2492 sage: J1 = HadamardEJA(2,QQ)
2493 sage: J2 = JordanSpinEJA(3,QQ)
2494 sage: J = DirectSumEJA(J1,J2)
2496 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2497 Euclidean Jordan algebra of dimension 3 over Rational Field)
2500 return self
._factors
2502 def projections(self
):
2504 Return a pair of projections onto this algebra's factors.
2508 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2509 ....: ComplexHermitianEJA,
2514 sage: J1 = JordanSpinEJA(2)
2515 sage: J2 = ComplexHermitianEJA(2)
2516 sage: J = DirectSumEJA(J1,J2)
2517 sage: (pi_left, pi_right) = J.projections()
2518 sage: J.one().to_vector()
2520 sage: pi_left(J.one()).to_vector()
2522 sage: pi_right(J.one()).to_vector()
2526 (J1
,J2
) = self
.factors()
2528 V_basis
= self
.vector_space().basis()
2529 P1
= matrix(self
.base_ring(), V_basis
[:n
])
2530 P2
= matrix(self
.base_ring(), V_basis
[n
:])
2531 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2532 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2533 return (pi_left
, pi_right
)
2535 def inclusions(self
):
2537 Return the pair of inclusion maps from our factors into us.
2541 sage: from mjo.eja.eja_algebra import (random_eja,
2542 ....: JordanSpinEJA,
2543 ....: RealSymmetricEJA,
2548 sage: J1 = JordanSpinEJA(3)
2549 sage: J2 = RealSymmetricEJA(2)
2550 sage: J = DirectSumEJA(J1,J2)
2551 sage: (iota_left, iota_right) = J.inclusions()
2552 sage: iota_left(J1.zero()) == J.zero()
2554 sage: iota_right(J2.zero()) == J.zero()
2556 sage: J1.one().to_vector()
2558 sage: iota_left(J1.one()).to_vector()
2560 sage: J2.one().to_vector()
2562 sage: iota_right(J2.one()).to_vector()
2564 sage: J.one().to_vector()
2569 Composing a projection with the corresponding inclusion should
2570 produce the identity map, and mismatching them should produce
2573 sage: set_random_seed()
2574 sage: J1 = random_eja()
2575 sage: J2 = random_eja()
2576 sage: J = DirectSumEJA(J1,J2)
2577 sage: (iota_left, iota_right) = J.inclusions()
2578 sage: (pi_left, pi_right) = J.projections()
2579 sage: pi_left*iota_left == J1.one().operator()
2581 sage: pi_right*iota_right == J2.one().operator()
2583 sage: (pi_left*iota_right).is_zero()
2585 sage: (pi_right*iota_left).is_zero()
2589 (J1
,J2
) = self
.factors()
2591 V_basis
= self
.vector_space().basis()
2592 I1
= matrix
.column(self
.base_ring(), V_basis
[:n
])
2593 I2
= matrix
.column(self
.base_ring(), V_basis
[n
:])
2594 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2595 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2596 return (iota_left
, iota_right
)
2598 def inner_product(self
, x
, y
):
2600 The standard Cartesian inner-product.
2602 We project ``x`` and ``y`` onto our factors, and add up the
2603 inner-products from the subalgebras.
2608 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2609 ....: QuaternionHermitianEJA,
2614 sage: J1 = HadamardEJA(3,QQ)
2615 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2616 sage: J = DirectSumEJA(J1,J2)
2621 sage: x1.inner_product(x2)
2623 sage: y1.inner_product(y2)
2625 sage: J.one().inner_product(J.one())
2629 (pi_left
, pi_right
) = self
.projections()
2635 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2639 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance