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_utils
import _mat2vec
38 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
40 def _coerce_map_from_base_ring(self
):
42 Disable the map from the base ring into the algebra.
44 Performing a nonsense conversion like this automatically
45 is counterpedagogical. The fallback is to try the usual
46 element constructor, which should also fail.
50 sage: from mjo.eja.eja_algebra import random_eja
54 sage: set_random_seed()
55 sage: J = random_eja()
57 Traceback (most recent call last):
59 ValueError: not a naturally-represented algebra element
75 sage: from mjo.eja.eja_algebra import (
76 ....: FiniteDimensionalEuclideanJordanAlgebra,
82 By definition, Jordan multiplication commutes::
84 sage: set_random_seed()
85 sage: J = random_eja()
86 sage: x,y = J.random_elements(2)
92 The ``field`` we're given must be real with ``check_field=True``::
94 sage: JordanSpinEJA(2,QQbar)
95 Traceback (most recent call last):
97 ValueError: scalar field is not real
99 The multiplication table must be square with ``check_axioms=True``::
101 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
102 Traceback (most recent call last):
104 ValueError: multiplication table is not square
108 if not field
.is_subring(RR
):
109 # Note: this does return true for the real algebraic
110 # field, the rationals, and any quadratic field where
111 # we've specified a real embedding.
112 raise ValueError("scalar field is not real")
114 # The multiplication table had better be square
117 if not all( len(l
) == n
for l
in mult_table
):
118 raise ValueError("multiplication table is not square")
120 self
._natural
_basis
= natural_basis
123 category
= MagmaticAlgebras(field
).FiniteDimensional()
124 category
= category
.WithBasis().Unital()
126 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
131 self
.print_options(bracket
='')
133 # The multiplication table we're given is necessarily in terms
134 # of vectors, because we don't have an algebra yet for
135 # anything to be an element of. However, it's faster in the
136 # long run to have the multiplication table be in terms of
137 # algebra elements. We do this after calling the superclass
138 # constructor so that from_vector() knows what to do.
139 self
._multiplication
_table
= [
140 list(map(lambda x
: self
.from_vector(x
), ls
))
145 if not self
._is
_commutative
():
146 raise ValueError("algebra is not commutative")
147 if not self
._is
_jordanian
():
148 raise ValueError("Jordan identity does not hold")
149 if not self
._inner
_product
_is
_associative
():
150 raise ValueError("inner product is not associative")
152 def _element_constructor_(self
, elt
):
154 Construct an element of this algebra from its natural
157 This gets called only after the parent element _call_ method
158 fails to find a coercion for the argument.
162 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
164 ....: RealSymmetricEJA)
168 The identity in `S^n` is converted to the identity in the EJA::
170 sage: J = RealSymmetricEJA(3)
171 sage: I = matrix.identity(QQ,3)
172 sage: J(I) == J.one()
175 This skew-symmetric matrix can't be represented in the EJA::
177 sage: J = RealSymmetricEJA(3)
178 sage: A = matrix(QQ,3, lambda i,j: i-j)
180 Traceback (most recent call last):
182 ArithmeticError: vector is not in free module
186 Ensure that we can convert any element of the two non-matrix
187 simple algebras (whose natural representations are their usual
188 vector representations) back and forth faithfully::
190 sage: set_random_seed()
191 sage: J = HadamardEJA.random_instance()
192 sage: x = J.random_element()
193 sage: J(x.to_vector().column()) == x
195 sage: J = JordanSpinEJA.random_instance()
196 sage: x = J.random_element()
197 sage: J(x.to_vector().column()) == x
201 msg
= "not a naturally-represented algebra element"
203 # The superclass implementation of random_element()
204 # needs to be able to coerce "0" into the algebra.
206 elif elt
in self
.base_ring():
207 # Ensure that no base ring -> algebra coercion is performed
208 # by this method. There's some stupidity in sage that would
209 # otherwise propagate to this method; for example, sage thinks
210 # that the integer 3 belongs to the space of 2-by-2 matrices.
211 raise ValueError(msg
)
213 natural_basis
= self
.natural_basis()
214 basis_space
= natural_basis
[0].matrix_space()
215 if elt
not in basis_space
:
216 raise ValueError(msg
)
218 # Thanks for nothing! Matrix spaces aren't vector spaces in
219 # Sage, so we have to figure out its natural-basis coordinates
220 # ourselves. We use the basis space's ring instead of the
221 # element's ring because the basis space might be an algebraic
222 # closure whereas the base ring of the 3-by-3 identity matrix
223 # could be QQ instead of QQbar.
224 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
225 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
226 coords
= W
.coordinate_vector(_mat2vec(elt
))
227 return self
.from_vector(coords
)
231 Return a string representation of ``self``.
235 sage: from mjo.eja.eja_algebra import JordanSpinEJA
239 Ensure that it says what we think it says::
241 sage: JordanSpinEJA(2, field=AA)
242 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
243 sage: JordanSpinEJA(3, field=RDF)
244 Euclidean Jordan algebra of dimension 3 over Real Double Field
247 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
248 return fmt
.format(self
.dimension(), self
.base_ring())
250 def product_on_basis(self
, i
, j
):
251 return self
._multiplication
_table
[i
][j
]
253 def _is_commutative(self
):
255 Whether or not this algebra's multiplication table is commutative.
257 This method should of course always return ``True``, unless
258 this algebra was constructed with ``check_axioms=False`` and
259 passed an invalid multiplication table.
261 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
262 for i
in range(self
.dimension())
263 for j
in range(self
.dimension()) )
265 def _is_jordanian(self
):
267 Whether or not this algebra's multiplication table respects the
268 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
270 We only check one arrangement of `x` and `y`, so for a
271 ``True`` result to be truly true, you should also check
272 :meth:`_is_commutative`. This method should of course always
273 return ``True``, unless this algebra was constructed with
274 ``check_axioms=False`` and passed an invalid multiplication table.
276 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
278 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
279 for i
in range(self
.dimension())
280 for j
in range(self
.dimension()) )
282 def _inner_product_is_associative(self
):
284 Return whether or not this algebra's inner product `B` is
285 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
287 This method should of course always return ``True``, unless
288 this algebra was constructed with ``check_axioms=False`` and
289 passed an invalid multiplication table.
292 # Used to check whether or not something is zero in an inexact
293 # ring. This number is sufficient to allow the construction of
294 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
297 for i
in range(self
.dimension()):
298 for j
in range(self
.dimension()):
299 for k
in range(self
.dimension()):
303 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
305 if self
.base_ring().is_exact():
309 if diff
.abs() > epsilon
:
315 def characteristic_polynomial_of(self
):
317 Return the algebra's "characteristic polynomial of" function,
318 which is itself a multivariate polynomial that, when evaluated
319 at the coordinates of some algebra element, returns that
320 element's characteristic polynomial.
322 The resulting polynomial has `n+1` variables, where `n` is the
323 dimension of this algebra. The first `n` variables correspond to
324 the coordinates of an algebra element: when evaluated at the
325 coordinates of an algebra element with respect to a certain
326 basis, the result is a univariate polynomial (in the one
327 remaining variable ``t``), namely the characteristic polynomial
332 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
336 The characteristic polynomial in the spin algebra is given in
337 Alizadeh, Example 11.11::
339 sage: J = JordanSpinEJA(3)
340 sage: p = J.characteristic_polynomial_of(); p
341 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
342 sage: xvec = J.one().to_vector()
346 By definition, the characteristic polynomial is a monic
347 degree-zero polynomial in a rank-zero algebra. Note that
348 Cayley-Hamilton is indeed satisfied since the polynomial
349 ``1`` evaluates to the identity element of the algebra on
352 sage: J = TrivialEJA()
353 sage: J.characteristic_polynomial_of()
360 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
361 a
= self
._charpoly
_coefficients
()
363 # We go to a bit of trouble here to reorder the
364 # indeterminates, so that it's easier to evaluate the
365 # characteristic polynomial at x's coordinates and get back
366 # something in terms of t, which is what we want.
367 S
= PolynomialRing(self
.base_ring(),'t')
371 S
= PolynomialRing(S
, R
.variable_names())
374 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
377 def inner_product(self
, x
, y
):
379 The inner product associated with this Euclidean Jordan algebra.
381 Defaults to the trace inner product, but can be overridden by
382 subclasses if they are sure that the necessary properties are
387 sage: from mjo.eja.eja_algebra import random_eja
391 Our inner product is "associative," which means the following for
392 a symmetric bilinear form::
394 sage: set_random_seed()
395 sage: J = random_eja()
396 sage: x,y,z = J.random_elements(3)
397 sage: (x*y).inner_product(z) == y.inner_product(x*z)
401 X
= x
.natural_representation()
402 Y
= y
.natural_representation()
403 return self
.natural_inner_product(X
,Y
)
406 def is_trivial(self
):
408 Return whether or not this algebra is trivial.
410 A trivial algebra contains only the zero element.
414 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
419 sage: J = ComplexHermitianEJA(3)
425 sage: J = TrivialEJA()
430 return self
.dimension() == 0
433 def multiplication_table(self
):
435 Return a visual representation of this algebra's multiplication
436 table (on basis elements).
440 sage: from mjo.eja.eja_algebra import JordanSpinEJA
444 sage: J = JordanSpinEJA(4)
445 sage: J.multiplication_table()
446 +----++----+----+----+----+
447 | * || e0 | e1 | e2 | e3 |
448 +====++====+====+====+====+
449 | e0 || e0 | e1 | e2 | e3 |
450 +----++----+----+----+----+
451 | e1 || e1 | e0 | 0 | 0 |
452 +----++----+----+----+----+
453 | e2 || e2 | 0 | e0 | 0 |
454 +----++----+----+----+----+
455 | e3 || e3 | 0 | 0 | e0 |
456 +----++----+----+----+----+
459 M
= list(self
._multiplication
_table
) # copy
460 for i
in range(len(M
)):
461 # M had better be "square"
462 M
[i
] = [self
.monomial(i
)] + M
[i
]
463 M
= [["*"] + list(self
.gens())] + M
464 return table(M
, header_row
=True, header_column
=True, frame
=True)
467 def natural_basis(self
):
469 Return a more-natural representation of this algebra's basis.
471 Every finite-dimensional Euclidean Jordan Algebra is a direct
472 sum of five simple algebras, four of which comprise Hermitian
473 matrices. This method returns the original "natural" basis
474 for our underlying vector space. (Typically, the natural basis
475 is used to construct the multiplication table in the first place.)
477 Note that this will always return a matrix. The standard basis
478 in `R^n` will be returned as `n`-by-`1` column matrices.
482 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
483 ....: RealSymmetricEJA)
487 sage: J = RealSymmetricEJA(2)
489 Finite family {0: e0, 1: e1, 2: e2}
490 sage: J.natural_basis()
492 [1 0] [ 0 0.7071067811865475?] [0 0]
493 [0 0], [0.7071067811865475? 0], [0 1]
498 sage: J = JordanSpinEJA(2)
500 Finite family {0: e0, 1: e1}
501 sage: J.natural_basis()
508 if self
._natural
_basis
is None:
509 M
= self
.natural_basis_space()
510 return tuple( M(b
.to_vector()) for b
in self
.basis() )
512 return self
._natural
_basis
515 def natural_basis_space(self
):
517 Return the matrix space in which this algebra's natural basis
520 Generally this will be an `n`-by-`1` column-vector space,
521 except when the algebra is trivial. There it's `n`-by-`n`
522 (where `n` is zero), to ensure that two elements of the
523 natural basis space (empty matrices) can be multiplied.
525 if self
.is_trivial():
526 return MatrixSpace(self
.base_ring(), 0)
527 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
528 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
530 return self
._natural
_basis
[0].matrix_space()
534 def natural_inner_product(X
,Y
):
536 Compute the inner product of two naturally-represented elements.
538 For example in the real symmetric matrix EJA, this will compute
539 the trace inner-product of two n-by-n symmetric matrices. The
540 default should work for the real cartesian product EJA, the
541 Jordan spin EJA, and the real symmetric matrices. The others
542 will have to be overridden.
544 return (X
.conjugate_transpose()*Y
).trace()
550 Return the unit element of this algebra.
554 sage: from mjo.eja.eja_algebra import (HadamardEJA,
559 sage: J = HadamardEJA(5)
561 e0 + e1 + e2 + e3 + e4
565 The identity element acts like the identity::
567 sage: set_random_seed()
568 sage: J = random_eja()
569 sage: x = J.random_element()
570 sage: J.one()*x == x and x*J.one() == x
573 The matrix of the unit element's operator is the identity::
575 sage: set_random_seed()
576 sage: J = random_eja()
577 sage: actual = J.one().operator().matrix()
578 sage: expected = matrix.identity(J.base_ring(), J.dimension())
579 sage: actual == expected
582 Ensure that the cached unit element (often precomputed by
583 hand) agrees with the computed one::
585 sage: set_random_seed()
586 sage: J = random_eja()
587 sage: cached = J.one()
588 sage: J.one.clear_cache()
589 sage: J.one() == cached
593 # We can brute-force compute the matrices of the operators
594 # that correspond to the basis elements of this algebra.
595 # If some linear combination of those basis elements is the
596 # algebra identity, then the same linear combination of
597 # their matrices has to be the identity matrix.
599 # Of course, matrices aren't vectors in sage, so we have to
600 # appeal to the "long vectors" isometry.
601 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
603 # Now we use basic linear algebra to find the coefficients,
604 # of the matrices-as-vectors-linear-combination, which should
605 # work for the original algebra basis too.
606 A
= matrix(self
.base_ring(), oper_vecs
)
608 # We used the isometry on the left-hand side already, but we
609 # still need to do it for the right-hand side. Recall that we
610 # wanted something that summed to the identity matrix.
611 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
613 # Now if there's an identity element in the algebra, this
614 # should work. We solve on the left to avoid having to
615 # transpose the matrix "A".
616 return self
.from_vector(A
.solve_left(b
))
619 def peirce_decomposition(self
, c
):
621 The Peirce decomposition of this algebra relative to the
624 In the future, this can be extended to a complete system of
625 orthogonal idempotents.
629 - ``c`` -- an idempotent of this algebra.
633 A triple (J0, J5, J1) containing two subalgebras and one subspace
636 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
637 corresponding to the eigenvalue zero.
639 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
640 corresponding to the eigenvalue one-half.
642 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
643 corresponding to the eigenvalue one.
645 These are the only possible eigenspaces for that operator, and this
646 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
647 orthogonal, and are subalgebras of this algebra with the appropriate
652 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
656 The canonical example comes from the symmetric matrices, which
657 decompose into diagonal and off-diagonal parts::
659 sage: J = RealSymmetricEJA(3)
660 sage: C = matrix(QQ, [ [1,0,0],
664 sage: J0,J5,J1 = J.peirce_decomposition(c)
666 Euclidean Jordan algebra of dimension 1...
668 Vector space of degree 6 and dimension 2...
670 Euclidean Jordan algebra of dimension 3...
671 sage: J0.one().natural_representation()
675 sage: orig_df = AA.options.display_format
676 sage: AA.options.display_format = 'radical'
677 sage: J.from_vector(J5.basis()[0]).natural_representation()
681 sage: J.from_vector(J5.basis()[1]).natural_representation()
685 sage: AA.options.display_format = orig_df
686 sage: J1.one().natural_representation()
693 Every algebra decomposes trivially with respect to its identity
696 sage: set_random_seed()
697 sage: J = random_eja()
698 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
699 sage: J0.dimension() == 0 and J5.dimension() == 0
701 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
704 The decomposition is into eigenspaces, and its components are
705 therefore necessarily orthogonal. Moreover, the identity
706 elements in the two subalgebras are the projections onto their
707 respective subspaces of the superalgebra's identity element::
709 sage: set_random_seed()
710 sage: J = random_eja()
711 sage: x = J.random_element()
712 sage: if not J.is_trivial():
713 ....: while x.is_nilpotent():
714 ....: x = J.random_element()
715 sage: c = x.subalgebra_idempotent()
716 sage: J0,J5,J1 = J.peirce_decomposition(c)
718 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
719 ....: w = w.superalgebra_element()
720 ....: y = J.from_vector(y)
721 ....: z = z.superalgebra_element()
722 ....: ipsum += w.inner_product(y).abs()
723 ....: ipsum += w.inner_product(z).abs()
724 ....: ipsum += y.inner_product(z).abs()
727 sage: J1(c) == J1.one()
729 sage: J0(J.one() - c) == J0.one()
733 if not c
.is_idempotent():
734 raise ValueError("element is not idempotent: %s" % c
)
736 # Default these to what they should be if they turn out to be
737 # trivial, because eigenspaces_left() won't return eigenvalues
738 # corresponding to trivial spaces (e.g. it returns only the
739 # eigenspace corresponding to lambda=1 if you take the
740 # decomposition relative to the identity element).
741 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
742 J0
= trivial
# eigenvalue zero
743 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
744 J1
= trivial
# eigenvalue one
746 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
747 if eigval
== ~
(self
.base_ring()(2)):
750 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
751 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
759 raise ValueError("unexpected eigenvalue: %s" % eigval
)
764 def random_element(self
, thorough
=False):
766 Return a random element of this algebra.
768 Our algebra superclass method only returns a linear
769 combination of at most two basis elements. We instead
770 want the vector space "random element" method that
771 returns a more diverse selection.
775 - ``thorough`` -- (boolean; default False) whether or not we
776 should generate irrational coefficients for the random
777 element when our base ring is irrational; this slows the
778 algebra operations to a crawl, but any truly random method
782 # For a general base ring... maybe we can trust this to do the
783 # right thing? Unlikely, but.
784 V
= self
.vector_space()
785 v
= V
.random_element()
787 if self
.base_ring() is AA
:
788 # The "random element" method of the algebraic reals is
789 # stupid at the moment, and only returns integers between
790 # -2 and 2, inclusive:
792 # https://trac.sagemath.org/ticket/30875
794 # Instead, we implement our own "random vector" method,
795 # and then coerce that into the algebra. We use the vector
796 # space degree here instead of the dimension because a
797 # subalgebra could (for example) be spanned by only two
798 # vectors, each with five coordinates. We need to
799 # generate all five coordinates.
801 v
*= QQbar
.random_element().real()
803 v
*= QQ
.random_element()
805 return self
.from_vector(V
.coordinate_vector(v
))
807 def random_elements(self
, count
, thorough
=False):
809 Return ``count`` random elements as a tuple.
813 - ``thorough`` -- (boolean; default False) whether or not we
814 should generate irrational coefficients for the random
815 elements when our base ring is irrational; this slows the
816 algebra operations to a crawl, but any truly random method
821 sage: from mjo.eja.eja_algebra import JordanSpinEJA
825 sage: J = JordanSpinEJA(3)
826 sage: x,y,z = J.random_elements(3)
827 sage: all( [ x in J, y in J, z in J ])
829 sage: len( J.random_elements(10) ) == 10
833 return tuple( self
.random_element(thorough
)
834 for idx
in range(count
) )
838 def _charpoly_coefficients(self
):
840 The `r` polynomial coefficients of the "characteristic polynomial
844 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
845 R
= PolynomialRing(self
.base_ring(), var_names
)
847 F
= R
.fraction_field()
850 # From a result in my book, these are the entries of the
851 # basis representation of L_x.
852 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
855 L_x
= matrix(F
, n
, n
, L_x_i_j
)
858 if self
.rank
.is_in_cache():
860 # There's no need to pad the system with redundant
861 # columns if we *know* they'll be redundant.
864 # Compute an extra power in case the rank is equal to
865 # the dimension (otherwise, we would stop at x^(r-1)).
866 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
867 for k
in range(n
+1) ]
868 A
= matrix
.column(F
, x_powers
[:n
])
869 AE
= A
.extended_echelon_form()
876 # The theory says that only the first "r" coefficients are
877 # nonzero, and they actually live in the original polynomial
878 # ring and not the fraction field. We negate them because
879 # in the actual characteristic polynomial, they get moved
880 # to the other side where x^r lives.
881 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
886 Return the rank of this EJA.
888 This is a cached method because we know the rank a priori for
889 all of the algebras we can construct. Thus we can avoid the
890 expensive ``_charpoly_coefficients()`` call unless we truly
891 need to compute the whole characteristic polynomial.
895 sage: from mjo.eja.eja_algebra import (HadamardEJA,
897 ....: RealSymmetricEJA,
898 ....: ComplexHermitianEJA,
899 ....: QuaternionHermitianEJA,
904 The rank of the Jordan spin algebra is always two::
906 sage: JordanSpinEJA(2).rank()
908 sage: JordanSpinEJA(3).rank()
910 sage: JordanSpinEJA(4).rank()
913 The rank of the `n`-by-`n` Hermitian real, complex, or
914 quaternion matrices is `n`::
916 sage: RealSymmetricEJA(4).rank()
918 sage: ComplexHermitianEJA(3).rank()
920 sage: QuaternionHermitianEJA(2).rank()
925 Ensure that every EJA that we know how to construct has a
926 positive integer rank, unless the algebra is trivial in
927 which case its rank will be zero::
929 sage: set_random_seed()
930 sage: J = random_eja()
934 sage: r > 0 or (r == 0 and J.is_trivial())
937 Ensure that computing the rank actually works, since the ranks
938 of all simple algebras are known and will be cached by default::
940 sage: J = HadamardEJA(4)
941 sage: J.rank.clear_cache()
947 sage: J = JordanSpinEJA(4)
948 sage: J.rank.clear_cache()
954 sage: J = RealSymmetricEJA(3)
955 sage: J.rank.clear_cache()
961 sage: J = ComplexHermitianEJA(2)
962 sage: J.rank.clear_cache()
968 sage: J = QuaternionHermitianEJA(2)
969 sage: J.rank.clear_cache()
973 return len(self
._charpoly
_coefficients
())
976 def vector_space(self
):
978 Return the vector space that underlies this algebra.
982 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
986 sage: J = RealSymmetricEJA(2)
987 sage: J.vector_space()
988 Vector space of dimension 3 over...
991 return self
.zero().to_vector().parent().ambient_vector_space()
994 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
997 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
999 Algebras whose basis consists of vectors with rational
1000 entries. Equivalently, algebras whose multiplication tables
1001 contain only rational coefficients.
1003 When an EJA has a basis that can be made rational, we can speed up
1004 the computation of its characteristic polynomial by doing it over
1005 ``QQ``. All of the named EJA constructors that we provide fall
1009 def _charpoly_coefficients(self
):
1011 Override the parent method with something that tries to compute
1012 over a faster (non-extension) field.
1016 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1020 The base ring of the resulting polynomial coefficients is what
1021 it should be, and not the rationals (unless the algebra was
1022 already over the rationals)::
1024 sage: J = JordanSpinEJA(3)
1025 sage: J._charpoly_coefficients()
1026 (X1^2 - X2^2 - X3^2, -2*X1)
1027 sage: a0 = J._charpoly_coefficients()[0]
1029 Algebraic Real Field
1030 sage: a0.base_ring()
1031 Algebraic Real Field
1034 if self
.base_ring() is QQ
:
1035 # There's no need to construct *another* algebra over the
1036 # rationals if this one is already over the rationals.
1037 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1038 return superclass
._charpoly
_coefficients
()
1041 map(lambda x
: x
.to_vector(), ls
)
1042 for ls
in self
._multiplication
_table
1045 # Do the computation over the rationals. The answer will be
1046 # the same, because our basis coordinates are (essentially)
1048 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1052 a
= J
._charpoly
_coefficients
()
1053 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1056 class ConcreteEuclideanJordanAlgebra
:
1058 A class for the Euclidean Jordan algebras that we know by name.
1060 These are the Jordan algebras whose basis, multiplication table,
1061 rank, and so on are known a priori. More to the point, they are
1062 the Euclidean Jordan algebras for which we are able to conjure up
1063 a "random instance."
1067 def _max_random_instance_size():
1069 Return an integer "size" that is an upper bound on the size of
1070 this algebra when it is used in a random test
1071 case. Unfortunately, the term "size" is ambiguous -- when
1072 dealing with `R^n` under either the Hadamard or Jordan spin
1073 product, the "size" refers to the dimension `n`. When dealing
1074 with a matrix algebra (real symmetric or complex/quaternion
1075 Hermitian), it refers to the size of the matrix, which is far
1076 less than the dimension of the underlying vector space.
1078 This method must be implemented in each subclass.
1080 raise NotImplementedError
1083 def random_instance(cls
, field
=AA
, **kwargs
):
1085 Return a random instance of this type of algebra.
1087 This method should be implemented in each subclass.
1089 from sage
.misc
.prandom
import choice
1090 eja_class
= choice(cls
.__subclasses
__())
1091 return eja_class
.random_instance(field
)
1094 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1096 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1098 Compared to the superclass constructor, we take a basis instead of
1099 a multiplication table because the latter can be computed in terms
1100 of the former when the product is known (like it is here).
1102 # Used in this class's fast _charpoly_coefficients() override.
1103 self
._basis
_normalizers
= None
1105 # We're going to loop through this a few times, so now's a good
1106 # time to ensure that it isn't a generator expression.
1107 basis
= tuple(basis
)
1109 algebra_dim
= len(basis
)
1110 if algebra_dim
> 1 and normalize_basis
:
1111 # We'll need sqrt(2) to normalize the basis, and this
1112 # winds up in the multiplication table, so the whole
1113 # algebra needs to be over the field extension.
1114 R
= PolynomialRing(field
, 'z')
1117 if p
.is_irreducible():
1118 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1119 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1120 self
._basis
_normalizers
= tuple(
1121 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1122 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1124 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1126 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1128 natural_basis
=basis
,
1131 if algebra_dim
== 0:
1132 self
.one
.set_cache(self
.zero())
1134 n
= basis
[0].nrows()
1135 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1136 # details of this algebra.
1137 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1141 def _charpoly_coefficients(self
):
1143 Override the parent method with something that tries to compute
1144 over a faster (non-extension) field.
1146 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1147 # We didn't normalize, or the basis we started with had
1148 # entries in a nice field already. Just compute the thing.
1149 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1151 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1152 self
._basis
_normalizers
) )
1154 # Do this over the rationals and convert back at the end.
1155 # Only works because we know the entries of the basis are
1156 # integers. The argument ``check_axioms=False`` is required
1157 # because the trace inner-product method for this
1158 # class is a stub and can't actually be checked.
1159 J
= MatrixEuclideanJordanAlgebra(QQ
,
1161 normalize_basis
=False,
1164 a
= J
._charpoly
_coefficients
()
1166 # Unfortunately, changing the basis does change the
1167 # coefficients of the characteristic polynomial, but since
1168 # these are really the coefficients of the "characteristic
1169 # polynomial of" function, everything is still nice and
1170 # unevaluated. It's therefore "obvious" how scaling the
1171 # basis affects the coordinate variables X1, X2, et
1172 # cetera. Scaling the first basis vector up by "n" adds a
1173 # factor of 1/n into every "X1" term, for example. So here
1174 # we simply undo the basis_normalizer scaling that we
1175 # performed earlier.
1177 # The a[0] access here is safe because trivial algebras
1178 # won't have any basis normalizers and therefore won't
1179 # make it to this "else" branch.
1180 XS
= a
[0].parent().gens()
1181 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1182 for i
in range(len(XS
)) }
1183 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1187 def multiplication_table_from_matrix_basis(basis
):
1189 At least three of the five simple Euclidean Jordan algebras have the
1190 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1191 multiplication on the right is matrix multiplication. Given a basis
1192 for the underlying matrix space, this function returns a
1193 multiplication table (obtained by looping through the basis
1194 elements) for an algebra of those matrices.
1196 # In S^2, for example, we nominally have four coordinates even
1197 # though the space is of dimension three only. The vector space V
1198 # is supposed to hold the entire long vector, and the subspace W
1199 # of V will be spanned by the vectors that arise from symmetric
1200 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1204 field
= basis
[0].base_ring()
1205 dimension
= basis
[0].nrows()
1207 V
= VectorSpace(field
, dimension
**2)
1208 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1210 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1213 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1214 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1222 Embed the matrix ``M`` into a space of real matrices.
1224 The matrix ``M`` can have entries in any field at the moment:
1225 the real numbers, complex numbers, or quaternions. And although
1226 they are not a field, we can probably support octonions at some
1227 point, too. This function returns a real matrix that "acts like"
1228 the original with respect to matrix multiplication; i.e.
1230 real_embed(M*N) = real_embed(M)*real_embed(N)
1233 raise NotImplementedError
1237 def real_unembed(M
):
1239 The inverse of :meth:`real_embed`.
1241 raise NotImplementedError
1245 def natural_inner_product(cls
,X
,Y
):
1246 Xu
= cls
.real_unembed(X
)
1247 Yu
= cls
.real_unembed(Y
)
1248 tr
= (Xu
*Yu
).trace()
1251 # Works in QQ, AA, RDF, et cetera.
1253 except AttributeError:
1254 # A quaternion doesn't have a real() method, but does
1255 # have coefficient_tuple() method that returns the
1256 # coefficients of 1, i, j, and k -- in that order.
1257 return tr
.coefficient_tuple()[0]
1260 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1264 The identity function, for embedding real matrices into real
1270 def real_unembed(M
):
1272 The identity function, for unembedding real matrices from real
1278 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1279 ConcreteEuclideanJordanAlgebra
):
1281 The rank-n simple EJA consisting of real symmetric n-by-n
1282 matrices, the usual symmetric Jordan product, and the trace inner
1283 product. It has dimension `(n^2 + n)/2` over the reals.
1287 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1291 sage: J = RealSymmetricEJA(2)
1292 sage: e0, e1, e2 = J.gens()
1300 In theory, our "field" can be any subfield of the reals::
1302 sage: RealSymmetricEJA(2, RDF)
1303 Euclidean Jordan algebra of dimension 3 over Real Double Field
1304 sage: RealSymmetricEJA(2, RR)
1305 Euclidean Jordan algebra of dimension 3 over Real Field with
1306 53 bits of precision
1310 The dimension of this algebra is `(n^2 + n) / 2`::
1312 sage: set_random_seed()
1313 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1314 sage: n = ZZ.random_element(1, n_max)
1315 sage: J = RealSymmetricEJA(n)
1316 sage: J.dimension() == (n^2 + n)/2
1319 The Jordan multiplication is what we think it is::
1321 sage: set_random_seed()
1322 sage: J = RealSymmetricEJA.random_instance()
1323 sage: x,y = J.random_elements(2)
1324 sage: actual = (x*y).natural_representation()
1325 sage: X = x.natural_representation()
1326 sage: Y = y.natural_representation()
1327 sage: expected = (X*Y + Y*X)/2
1328 sage: actual == expected
1330 sage: J(expected) == x*y
1333 We can change the generator prefix::
1335 sage: RealSymmetricEJA(3, prefix='q').gens()
1336 (q0, q1, q2, q3, q4, q5)
1338 Our natural basis is normalized with respect to the natural inner
1339 product unless we specify otherwise::
1341 sage: set_random_seed()
1342 sage: J = RealSymmetricEJA.random_instance()
1343 sage: all( b.norm() == 1 for b in J.gens() )
1346 Since our natural basis is normalized with respect to the natural
1347 inner product, and since we know that this algebra is an EJA, any
1348 left-multiplication operator's matrix will be symmetric because
1349 natural->EJA basis representation is an isometry and within the EJA
1350 the operator is self-adjoint by the Jordan axiom::
1352 sage: set_random_seed()
1353 sage: x = RealSymmetricEJA.random_instance().random_element()
1354 sage: x.operator().matrix().is_symmetric()
1357 We can construct the (trivial) algebra of rank zero::
1359 sage: RealSymmetricEJA(0)
1360 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1364 def _denormalized_basis(cls
, n
, field
):
1366 Return a basis for the space of real symmetric n-by-n matrices.
1370 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1374 sage: set_random_seed()
1375 sage: n = ZZ.random_element(1,5)
1376 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1377 sage: all( M.is_symmetric() for M in B)
1381 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1385 for j
in range(i
+1):
1386 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1390 Sij
= Eij
+ Eij
.transpose()
1396 def _max_random_instance_size():
1397 return 4 # Dimension 10
1400 def random_instance(cls
, field
=AA
, **kwargs
):
1402 Return a random instance of this type of algebra.
1404 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1405 return cls(n
, field
, **kwargs
)
1407 def __init__(self
, n
, field
=AA
, **kwargs
):
1408 basis
= self
._denormalized
_basis
(n
, field
)
1409 super(RealSymmetricEJA
, self
).__init
__(field
,
1413 self
.rank
.set_cache(n
)
1416 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1420 Embed the n-by-n complex matrix ``M`` into the space of real
1421 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1422 bi` to the block matrix ``[[a,b],[-b,a]]``.
1426 sage: from mjo.eja.eja_algebra import \
1427 ....: ComplexMatrixEuclideanJordanAlgebra
1431 sage: F = QuadraticField(-1, 'I')
1432 sage: x1 = F(4 - 2*i)
1433 sage: x2 = F(1 + 2*i)
1436 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1437 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1446 Embedding is a homomorphism (isomorphism, in fact)::
1448 sage: set_random_seed()
1449 sage: n = ZZ.random_element(3)
1450 sage: F = QuadraticField(-1, 'I')
1451 sage: X = random_matrix(F, n)
1452 sage: Y = random_matrix(F, n)
1453 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1454 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1455 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1462 raise ValueError("the matrix 'M' must be square")
1464 # We don't need any adjoined elements...
1465 field
= M
.base_ring().base_ring()
1469 a
= z
.list()[0] # real part, I guess
1470 b
= z
.list()[1] # imag part, I guess
1471 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1473 return matrix
.block(field
, n
, blocks
)
1477 def real_unembed(M
):
1479 The inverse of _embed_complex_matrix().
1483 sage: from mjo.eja.eja_algebra import \
1484 ....: ComplexMatrixEuclideanJordanAlgebra
1488 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1489 ....: [-2, 1, -4, 3],
1490 ....: [ 9, 10, 11, 12],
1491 ....: [-10, 9, -12, 11] ])
1492 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1494 [ 10*I + 9 12*I + 11]
1498 Unembedding is the inverse of embedding::
1500 sage: set_random_seed()
1501 sage: F = QuadraticField(-1, 'I')
1502 sage: M = random_matrix(F, 3)
1503 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1504 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1510 raise ValueError("the matrix 'M' must be square")
1511 if not n
.mod(2).is_zero():
1512 raise ValueError("the matrix 'M' must be a complex embedding")
1514 # If "M" was normalized, its base ring might have roots
1515 # adjoined and they can stick around after unembedding.
1516 field
= M
.base_ring()
1517 R
= PolynomialRing(field
, 'z')
1520 # Sage doesn't know how to embed AA into QQbar, i.e. how
1521 # to adjoin sqrt(-1) to AA.
1524 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1527 # Go top-left to bottom-right (reading order), converting every
1528 # 2-by-2 block we see to a single complex element.
1530 for k
in range(n
/2):
1531 for j
in range(n
/2):
1532 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1533 if submat
[0,0] != submat
[1,1]:
1534 raise ValueError('bad on-diagonal submatrix')
1535 if submat
[0,1] != -submat
[1,0]:
1536 raise ValueError('bad off-diagonal submatrix')
1537 z
= submat
[0,0] + submat
[0,1]*i
1540 return matrix(F
, n
/2, elements
)
1544 def natural_inner_product(cls
,X
,Y
):
1546 Compute a natural inner product in this algebra directly from
1551 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1555 This gives the same answer as the slow, default method implemented
1556 in :class:`MatrixEuclideanJordanAlgebra`::
1558 sage: set_random_seed()
1559 sage: J = ComplexHermitianEJA.random_instance()
1560 sage: x,y = J.random_elements(2)
1561 sage: Xe = x.natural_representation()
1562 sage: Ye = y.natural_representation()
1563 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1564 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1565 sage: expected = (X*Y).trace().real()
1566 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1567 sage: actual == expected
1571 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1574 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1575 ConcreteEuclideanJordanAlgebra
):
1577 The rank-n simple EJA consisting of complex Hermitian n-by-n
1578 matrices over the real numbers, the usual symmetric Jordan product,
1579 and the real-part-of-trace inner product. It has dimension `n^2` over
1584 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1588 In theory, our "field" can be any subfield of the reals::
1590 sage: ComplexHermitianEJA(2, RDF)
1591 Euclidean Jordan algebra of dimension 4 over Real Double Field
1592 sage: ComplexHermitianEJA(2, RR)
1593 Euclidean Jordan algebra of dimension 4 over Real Field with
1594 53 bits of precision
1598 The dimension of this algebra is `n^2`::
1600 sage: set_random_seed()
1601 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1602 sage: n = ZZ.random_element(1, n_max)
1603 sage: J = ComplexHermitianEJA(n)
1604 sage: J.dimension() == n^2
1607 The Jordan multiplication is what we think it is::
1609 sage: set_random_seed()
1610 sage: J = ComplexHermitianEJA.random_instance()
1611 sage: x,y = J.random_elements(2)
1612 sage: actual = (x*y).natural_representation()
1613 sage: X = x.natural_representation()
1614 sage: Y = y.natural_representation()
1615 sage: expected = (X*Y + Y*X)/2
1616 sage: actual == expected
1618 sage: J(expected) == x*y
1621 We can change the generator prefix::
1623 sage: ComplexHermitianEJA(2, prefix='z').gens()
1626 Our natural basis is normalized with respect to the natural inner
1627 product unless we specify otherwise::
1629 sage: set_random_seed()
1630 sage: J = ComplexHermitianEJA.random_instance()
1631 sage: all( b.norm() == 1 for b in J.gens() )
1634 Since our natural basis is normalized with respect to the natural
1635 inner product, and since we know that this algebra is an EJA, any
1636 left-multiplication operator's matrix will be symmetric because
1637 natural->EJA basis representation is an isometry and within the EJA
1638 the operator is self-adjoint by the Jordan axiom::
1640 sage: set_random_seed()
1641 sage: x = ComplexHermitianEJA.random_instance().random_element()
1642 sage: x.operator().matrix().is_symmetric()
1645 We can construct the (trivial) algebra of rank zero::
1647 sage: ComplexHermitianEJA(0)
1648 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1653 def _denormalized_basis(cls
, n
, field
):
1655 Returns a basis for the space of complex Hermitian n-by-n matrices.
1657 Why do we embed these? Basically, because all of numerical linear
1658 algebra assumes that you're working with vectors consisting of `n`
1659 entries from a field and scalars from the same field. There's no way
1660 to tell SageMath that (for example) the vectors contain complex
1661 numbers, while the scalar field is real.
1665 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1669 sage: set_random_seed()
1670 sage: n = ZZ.random_element(1,5)
1671 sage: field = QuadraticField(2, 'sqrt2')
1672 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1673 sage: all( M.is_symmetric() for M in B)
1677 R
= PolynomialRing(field
, 'z')
1679 F
= field
.extension(z
**2 + 1, 'I')
1682 # This is like the symmetric case, but we need to be careful:
1684 # * We want conjugate-symmetry, not just symmetry.
1685 # * The diagonal will (as a result) be real.
1689 for j
in range(i
+1):
1690 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1692 Sij
= cls
.real_embed(Eij
)
1695 # The second one has a minus because it's conjugated.
1696 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1698 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1701 # Since we embedded these, we can drop back to the "field" that we
1702 # started with instead of the complex extension "F".
1703 return ( s
.change_ring(field
) for s
in S
)
1706 def __init__(self
, n
, field
=AA
, **kwargs
):
1707 basis
= self
._denormalized
_basis
(n
,field
)
1708 super(ComplexHermitianEJA
,self
).__init
__(field
,
1712 self
.rank
.set_cache(n
)
1715 def _max_random_instance_size():
1716 return 3 # Dimension 9
1719 def random_instance(cls
, field
=AA
, **kwargs
):
1721 Return a random instance of this type of algebra.
1723 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1724 return cls(n
, field
, **kwargs
)
1726 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1730 Embed the n-by-n quaternion matrix ``M`` into the space of real
1731 matrices of size 4n-by-4n by first sending each quaternion entry `z
1732 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1733 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1738 sage: from mjo.eja.eja_algebra import \
1739 ....: QuaternionMatrixEuclideanJordanAlgebra
1743 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1744 sage: i,j,k = Q.gens()
1745 sage: x = 1 + 2*i + 3*j + 4*k
1746 sage: M = matrix(Q, 1, [[x]])
1747 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1753 Embedding is a homomorphism (isomorphism, in fact)::
1755 sage: set_random_seed()
1756 sage: n = ZZ.random_element(2)
1757 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1758 sage: X = random_matrix(Q, n)
1759 sage: Y = random_matrix(Q, n)
1760 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1761 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1762 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1767 quaternions
= M
.base_ring()
1770 raise ValueError("the matrix 'M' must be square")
1772 F
= QuadraticField(-1, 'I')
1777 t
= z
.coefficient_tuple()
1782 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1783 [-c
+ d
*i
, a
- b
*i
]])
1784 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1785 blocks
.append(realM
)
1787 # We should have real entries by now, so use the realest field
1788 # we've got for the return value.
1789 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1794 def real_unembed(M
):
1796 The inverse of _embed_quaternion_matrix().
1800 sage: from mjo.eja.eja_algebra import \
1801 ....: QuaternionMatrixEuclideanJordanAlgebra
1805 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1806 ....: [-2, 1, -4, 3],
1807 ....: [-3, 4, 1, -2],
1808 ....: [-4, -3, 2, 1]])
1809 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1810 [1 + 2*i + 3*j + 4*k]
1814 Unembedding is the inverse of embedding::
1816 sage: set_random_seed()
1817 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1818 sage: M = random_matrix(Q, 3)
1819 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1820 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1826 raise ValueError("the matrix 'M' must be square")
1827 if not n
.mod(4).is_zero():
1828 raise ValueError("the matrix 'M' must be a quaternion embedding")
1830 # Use the base ring of the matrix to ensure that its entries can be
1831 # multiplied by elements of the quaternion algebra.
1832 field
= M
.base_ring()
1833 Q
= QuaternionAlgebra(field
,-1,-1)
1836 # Go top-left to bottom-right (reading order), converting every
1837 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1840 for l
in range(n
/4):
1841 for m
in range(n
/4):
1842 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1843 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1844 if submat
[0,0] != submat
[1,1].conjugate():
1845 raise ValueError('bad on-diagonal submatrix')
1846 if submat
[0,1] != -submat
[1,0].conjugate():
1847 raise ValueError('bad off-diagonal submatrix')
1848 z
= submat
[0,0].real()
1849 z
+= submat
[0,0].imag()*i
1850 z
+= submat
[0,1].real()*j
1851 z
+= submat
[0,1].imag()*k
1854 return matrix(Q
, n
/4, elements
)
1858 def natural_inner_product(cls
,X
,Y
):
1860 Compute a natural inner product in this algebra directly from
1865 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1869 This gives the same answer as the slow, default method implemented
1870 in :class:`MatrixEuclideanJordanAlgebra`::
1872 sage: set_random_seed()
1873 sage: J = QuaternionHermitianEJA.random_instance()
1874 sage: x,y = J.random_elements(2)
1875 sage: Xe = x.natural_representation()
1876 sage: Ye = y.natural_representation()
1877 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1878 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1879 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1880 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1881 sage: actual == expected
1885 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1888 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1889 ConcreteEuclideanJordanAlgebra
):
1891 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1892 matrices, the usual symmetric Jordan product, and the
1893 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1898 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1902 In theory, our "field" can be any subfield of the reals::
1904 sage: QuaternionHermitianEJA(2, RDF)
1905 Euclidean Jordan algebra of dimension 6 over Real Double Field
1906 sage: QuaternionHermitianEJA(2, RR)
1907 Euclidean Jordan algebra of dimension 6 over Real Field with
1908 53 bits of precision
1912 The dimension of this algebra is `2*n^2 - n`::
1914 sage: set_random_seed()
1915 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
1916 sage: n = ZZ.random_element(1, n_max)
1917 sage: J = QuaternionHermitianEJA(n)
1918 sage: J.dimension() == 2*(n^2) - n
1921 The Jordan multiplication is what we think it is::
1923 sage: set_random_seed()
1924 sage: J = QuaternionHermitianEJA.random_instance()
1925 sage: x,y = J.random_elements(2)
1926 sage: actual = (x*y).natural_representation()
1927 sage: X = x.natural_representation()
1928 sage: Y = y.natural_representation()
1929 sage: expected = (X*Y + Y*X)/2
1930 sage: actual == expected
1932 sage: J(expected) == x*y
1935 We can change the generator prefix::
1937 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1938 (a0, a1, a2, a3, a4, a5)
1940 Our natural basis is normalized with respect to the natural inner
1941 product unless we specify otherwise::
1943 sage: set_random_seed()
1944 sage: J = QuaternionHermitianEJA.random_instance()
1945 sage: all( b.norm() == 1 for b in J.gens() )
1948 Since our natural basis is normalized with respect to the natural
1949 inner product, and since we know that this algebra is an EJA, any
1950 left-multiplication operator's matrix will be symmetric because
1951 natural->EJA basis representation is an isometry and within the EJA
1952 the operator is self-adjoint by the Jordan axiom::
1954 sage: set_random_seed()
1955 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1956 sage: x.operator().matrix().is_symmetric()
1959 We can construct the (trivial) algebra of rank zero::
1961 sage: QuaternionHermitianEJA(0)
1962 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1966 def _denormalized_basis(cls
, n
, field
):
1968 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1970 Why do we embed these? Basically, because all of numerical
1971 linear algebra assumes that you're working with vectors consisting
1972 of `n` entries from a field and scalars from the same field. There's
1973 no way to tell SageMath that (for example) the vectors contain
1974 complex numbers, while the scalar field is real.
1978 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1982 sage: set_random_seed()
1983 sage: n = ZZ.random_element(1,5)
1984 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1985 sage: all( M.is_symmetric() for M in B )
1989 Q
= QuaternionAlgebra(QQ
,-1,-1)
1992 # This is like the symmetric case, but we need to be careful:
1994 # * We want conjugate-symmetry, not just symmetry.
1995 # * The diagonal will (as a result) be real.
1999 for j
in range(i
+1):
2000 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2002 Sij
= cls
.real_embed(Eij
)
2005 # The second, third, and fourth ones have a minus
2006 # because they're conjugated.
2007 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2009 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2011 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2013 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2016 # Since we embedded these, we can drop back to the "field" that we
2017 # started with instead of the quaternion algebra "Q".
2018 return ( s
.change_ring(field
) for s
in S
)
2021 def __init__(self
, n
, field
=AA
, **kwargs
):
2022 basis
= self
._denormalized
_basis
(n
,field
)
2023 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2027 self
.rank
.set_cache(n
)
2030 def _max_random_instance_size():
2032 The maximum rank of a random QuaternionHermitianEJA.
2034 return 2 # Dimension 6
2037 def random_instance(cls
, field
=AA
, **kwargs
):
2039 Return a random instance of this type of algebra.
2041 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2042 return cls(n
, field
, **kwargs
)
2045 class HadamardEJA(RationalBasisEuclideanJordanAlgebra
,
2046 ConcreteEuclideanJordanAlgebra
):
2048 Return the Euclidean Jordan Algebra corresponding to the set
2049 `R^n` under the Hadamard product.
2051 Note: this is nothing more than the Cartesian product of ``n``
2052 copies of the spin algebra. Once Cartesian product algebras
2053 are implemented, this can go.
2057 sage: from mjo.eja.eja_algebra import HadamardEJA
2061 This multiplication table can be verified by hand::
2063 sage: J = HadamardEJA(3)
2064 sage: e0,e1,e2 = J.gens()
2080 We can change the generator prefix::
2082 sage: HadamardEJA(3, prefix='r').gens()
2086 def __init__(self
, n
, field
=AA
, **kwargs
):
2087 V
= VectorSpace(field
, n
)
2088 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
2091 super(HadamardEJA
, self
).__init
__(field
,
2095 self
.rank
.set_cache(n
)
2098 self
.one
.set_cache( self
.zero() )
2100 self
.one
.set_cache( sum(self
.gens()) )
2103 def _max_random_instance_size():
2105 The maximum dimension of a random HadamardEJA.
2110 def random_instance(cls
, field
=AA
, **kwargs
):
2112 Return a random instance of this type of algebra.
2114 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2115 return cls(n
, field
, **kwargs
)
2118 def inner_product(self
, x
, y
):
2120 Faster to reimplement than to use natural representations.
2124 sage: from mjo.eja.eja_algebra import HadamardEJA
2128 Ensure that this is the usual inner product for the algebras
2131 sage: set_random_seed()
2132 sage: J = HadamardEJA.random_instance()
2133 sage: x,y = J.random_elements(2)
2134 sage: X = x.natural_representation()
2135 sage: Y = y.natural_representation()
2136 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2140 return x
.to_vector().inner_product(y
.to_vector())
2143 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
,
2144 ConcreteEuclideanJordanAlgebra
):
2146 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2147 with the half-trace inner product and jordan product ``x*y =
2148 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2149 a symmetric positive-definite "bilinear form" matrix. Its
2150 dimension is the size of `B`, and it has rank two in dimensions
2151 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2152 the identity matrix of order ``n``.
2154 We insist that the one-by-one upper-left identity block of `B` be
2155 passed in as well so that we can be passed a matrix of size zero
2156 to construct a trivial algebra.
2160 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2161 ....: JordanSpinEJA)
2165 When no bilinear form is specified, the identity matrix is used,
2166 and the resulting algebra is the Jordan spin algebra::
2168 sage: B = matrix.identity(AA,3)
2169 sage: J0 = BilinearFormEJA(B)
2170 sage: J1 = JordanSpinEJA(3)
2171 sage: J0.multiplication_table() == J0.multiplication_table()
2174 An error is raised if the matrix `B` does not correspond to a
2175 positive-definite bilinear form::
2177 sage: B = matrix.random(QQ,2,3)
2178 sage: J = BilinearFormEJA(B)
2179 Traceback (most recent call last):
2181 ValueError: bilinear form is not positive-definite
2182 sage: B = matrix.zero(QQ,3)
2183 sage: J = BilinearFormEJA(B)
2184 Traceback (most recent call last):
2186 ValueError: bilinear form is not positive-definite
2190 We can create a zero-dimensional algebra::
2192 sage: B = matrix.identity(AA,0)
2193 sage: J = BilinearFormEJA(B)
2197 We can check the multiplication condition given in the Jordan, von
2198 Neumann, and Wigner paper (and also discussed on my "On the
2199 symmetry..." paper). Note that this relies heavily on the standard
2200 choice of basis, as does anything utilizing the bilinear form matrix::
2202 sage: set_random_seed()
2203 sage: n = ZZ.random_element(5)
2204 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2205 sage: B11 = matrix.identity(QQ,1)
2206 sage: B22 = M.transpose()*M
2207 sage: B = block_matrix(2,2,[ [B11,0 ],
2209 sage: J = BilinearFormEJA(B)
2210 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2211 sage: V = J.vector_space()
2212 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2213 ....: for ei in eis ]
2214 sage: actual = [ sis[i]*sis[j]
2215 ....: for i in range(n-1)
2216 ....: for j in range(n-1) ]
2217 sage: expected = [ J.one() if i == j else J.zero()
2218 ....: for i in range(n-1)
2219 ....: for j in range(n-1) ]
2220 sage: actual == expected
2223 def __init__(self
, B
, field
=AA
, **kwargs
):
2227 if not B
.is_positive_definite():
2228 raise ValueError("bilinear form is not positive-definite")
2230 V
= VectorSpace(field
, n
)
2231 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2240 z0
= (B
*x
).inner_product(y
)
2241 zbar
= y0
*xbar
+ x0
*ybar
2242 z
= V([z0
] + zbar
.list())
2243 mult_table
[i
][j
] = z
2245 # The rank of this algebra is two, unless we're in a
2246 # one-dimensional ambient space (because the rank is bounded
2247 # by the ambient dimension).
2248 super(BilinearFormEJA
, self
).__init
__(field
,
2252 self
.rank
.set_cache(min(n
,2))
2255 self
.one
.set_cache( self
.zero() )
2257 self
.one
.set_cache( self
.monomial(0) )
2260 def _max_random_instance_size():
2262 The maximum dimension of a random BilinearFormEJA.
2267 def random_instance(cls
, field
=AA
, **kwargs
):
2269 Return a random instance of this algebra.
2271 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2273 B
= matrix
.identity(field
, n
)
2274 return cls(B
, field
, **kwargs
)
2276 B11
= matrix
.identity(field
,1)
2277 M
= matrix
.random(field
, n
-1)
2278 I
= matrix
.identity(field
, n
-1)
2279 alpha
= field
.zero()
2280 while alpha
.is_zero():
2281 alpha
= field
.random_element().abs()
2282 B22
= M
.transpose()*M
+ alpha
*I
2284 from sage
.matrix
.special
import block_matrix
2285 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2288 return cls(B
, field
, **kwargs
)
2290 def inner_product(self
, x
, y
):
2292 Half of the trace inner product.
2294 This is defined so that the special case of the Jordan spin
2295 algebra gets the usual inner product.
2299 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2303 Ensure that this is one-half of the trace inner-product when
2304 the algebra isn't just the reals (when ``n`` isn't one). This
2305 is in Faraut and Koranyi, and also my "On the symmetry..."
2308 sage: set_random_seed()
2309 sage: J = BilinearFormEJA.random_instance()
2310 sage: n = J.dimension()
2311 sage: x = J.random_element()
2312 sage: y = J.random_element()
2313 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
2317 return (self
._B
*x
.to_vector()).inner_product(y
.to_vector())
2320 class JordanSpinEJA(BilinearFormEJA
):
2322 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2323 with the usual inner product and jordan product ``x*y =
2324 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2329 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2333 This multiplication table can be verified by hand::
2335 sage: J = JordanSpinEJA(4)
2336 sage: e0,e1,e2,e3 = J.gens()
2352 We can change the generator prefix::
2354 sage: JordanSpinEJA(2, prefix='B').gens()
2359 Ensure that we have the usual inner product on `R^n`::
2361 sage: set_random_seed()
2362 sage: J = JordanSpinEJA.random_instance()
2363 sage: x,y = J.random_elements(2)
2364 sage: X = x.natural_representation()
2365 sage: Y = y.natural_representation()
2366 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2370 def __init__(self
, n
, field
=AA
, **kwargs
):
2371 # This is a special case of the BilinearFormEJA with the identity
2372 # matrix as its bilinear form.
2373 B
= matrix
.identity(field
, n
)
2374 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2377 def _max_random_instance_size():
2379 The maximum dimension of a random JordanSpinEJA.
2384 def random_instance(cls
, field
=AA
, **kwargs
):
2386 Return a random instance of this type of algebra.
2388 Needed here to override the implementation for ``BilinearFormEJA``.
2390 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2391 return cls(n
, field
, **kwargs
)
2394 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
,
2395 ConcreteEuclideanJordanAlgebra
):
2397 The trivial Euclidean Jordan algebra consisting of only a zero element.
2401 sage: from mjo.eja.eja_algebra import TrivialEJA
2405 sage: J = TrivialEJA()
2412 sage: 7*J.one()*12*J.one()
2414 sage: J.one().inner_product(J.one())
2416 sage: J.one().norm()
2418 sage: J.one().subalgebra_generated_by()
2419 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2424 def __init__(self
, field
=AA
, **kwargs
):
2426 super(TrivialEJA
, self
).__init
__(field
,
2430 # The rank is zero using my definition, namely the dimension of the
2431 # largest subalgebra generated by any element.
2432 self
.rank
.set_cache(0)
2433 self
.one
.set_cache( self
.zero() )
2436 def random_instance(cls
, field
=AA
, **kwargs
):
2437 # We don't take a "size" argument so the superclass method is
2438 # inappropriate for us.
2439 return cls(field
, **kwargs
)
2441 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2443 The external (orthogonal) direct sum of two other Euclidean Jordan
2444 algebras. Essentially the Cartesian product of its two factors.
2445 Every Euclidean Jordan algebra decomposes into an orthogonal
2446 direct sum of simple Euclidean Jordan algebras, so no generality
2447 is lost by providing only this construction.
2451 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2452 ....: RealSymmetricEJA,
2457 sage: J1 = HadamardEJA(2)
2458 sage: J2 = RealSymmetricEJA(3)
2459 sage: J = DirectSumEJA(J1,J2)
2466 def __init__(self
, J1
, J2
, field
=AA
, **kwargs
):
2467 self
._factors
= (J1
, J2
)
2471 V
= VectorSpace(field
, n
)
2472 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2476 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2477 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2481 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2482 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2484 super(DirectSumEJA
, self
).__init
__(field
,
2488 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2493 Return the pair of this algebra's factors.
2497 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2498 ....: JordanSpinEJA,
2503 sage: J1 = HadamardEJA(2,QQ)
2504 sage: J2 = JordanSpinEJA(3,QQ)
2505 sage: J = DirectSumEJA(J1,J2)
2507 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2508 Euclidean Jordan algebra of dimension 3 over Rational Field)
2511 return self
._factors
2513 def projections(self
):
2515 Return a pair of projections onto this algebra's factors.
2519 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2520 ....: ComplexHermitianEJA,
2525 sage: J1 = JordanSpinEJA(2)
2526 sage: J2 = ComplexHermitianEJA(2)
2527 sage: J = DirectSumEJA(J1,J2)
2528 sage: (pi_left, pi_right) = J.projections()
2529 sage: J.one().to_vector()
2531 sage: pi_left(J.one()).to_vector()
2533 sage: pi_right(J.one()).to_vector()
2537 (J1
,J2
) = self
.factors()
2539 pi_left
= lambda x
: J1
.from_vector(x
.to_vector()[:n
])
2540 pi_right
= lambda x
: J2
.from_vector(x
.to_vector()[n
:])
2541 return (pi_left
, pi_right
)
2543 def inclusions(self
):
2545 Return the pair of inclusion maps from our factors into us.
2549 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2550 ....: RealSymmetricEJA,
2555 sage: J1 = JordanSpinEJA(3)
2556 sage: J2 = RealSymmetricEJA(2)
2557 sage: J = DirectSumEJA(J1,J2)
2558 sage: (iota_left, iota_right) = J.inclusions()
2559 sage: iota_left(J1.zero()) == J.zero()
2561 sage: iota_right(J2.zero()) == J.zero()
2563 sage: J1.one().to_vector()
2565 sage: iota_left(J1.one()).to_vector()
2567 sage: J2.one().to_vector()
2569 sage: iota_right(J2.one()).to_vector()
2571 sage: J.one().to_vector()
2575 (J1
,J2
) = self
.factors()
2577 V_basis
= self
.vector_space().basis()
2578 I1
= matrix
.column(self
.base_ring(), V_basis
[:n
])
2579 I2
= matrix
.column(self
.base_ring(), V_basis
[n
:])
2580 iota_left
= lambda x
: self
.from_vector(I1
*x
.to_vector())
2581 iota_right
= lambda x
: self
.from_vector(I2
*+x
.to_vector())
2582 return (iota_left
, iota_right
)
2584 def inner_product(self
, x
, y
):
2586 The standard Cartesian inner-product.
2588 We project ``x`` and ``y`` onto our factors, and add up the
2589 inner-products from the subalgebras.
2594 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2595 ....: QuaternionHermitianEJA,
2600 sage: J1 = HadamardEJA(3)
2601 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2602 sage: J = DirectSumEJA(J1,J2)
2607 sage: x1.inner_product(x2)
2609 sage: y1.inner_product(y2)
2611 sage: J.one().inner_product(J.one())
2615 (pi_left
, pi_right
) = self
.projections()
2621 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2625 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance