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,
390 ....: BilinearFormEJA)
394 Our inner product is "associative," which means the following for
395 a symmetric bilinear form::
397 sage: set_random_seed()
398 sage: J = random_eja()
399 sage: x,y,z = J.random_elements(3)
400 sage: (x*y).inner_product(z) == y.inner_product(x*z)
405 Ensure that this is the usual inner product for the algebras
408 sage: set_random_seed()
409 sage: J = HadamardEJA.random_instance()
410 sage: x,y = J.random_elements(2)
411 sage: actual = x.inner_product(y)
412 sage: expected = x.to_vector().inner_product(y.to_vector())
413 sage: actual == expected
416 Ensure that this is one-half of the trace inner-product in a
417 BilinearFormEJA that isn't just the reals (when ``n`` isn't
418 one). This is in Faraut and Koranyi, and also my "On the
421 sage: set_random_seed()
422 sage: J = BilinearFormEJA.random_instance()
423 sage: n = J.dimension()
424 sage: x = J.random_element()
425 sage: y = J.random_element()
426 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
429 B
= self
._inner
_product
_matrix
430 return (B
*x
.to_vector()).inner_product(y
.to_vector())
433 def is_trivial(self
):
435 Return whether or not this algebra is trivial.
437 A trivial algebra contains only the zero element.
441 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
446 sage: J = ComplexHermitianEJA(3)
452 sage: J = TrivialEJA()
457 return self
.dimension() == 0
460 def multiplication_table(self
):
462 Return a visual representation of this algebra's multiplication
463 table (on basis elements).
467 sage: from mjo.eja.eja_algebra import JordanSpinEJA
471 sage: J = JordanSpinEJA(4)
472 sage: J.multiplication_table()
473 +----++----+----+----+----+
474 | * || e0 | e1 | e2 | e3 |
475 +====++====+====+====+====+
476 | e0 || e0 | e1 | e2 | e3 |
477 +----++----+----+----+----+
478 | e1 || e1 | e0 | 0 | 0 |
479 +----++----+----+----+----+
480 | e2 || e2 | 0 | e0 | 0 |
481 +----++----+----+----+----+
482 | e3 || e3 | 0 | 0 | e0 |
483 +----++----+----+----+----+
486 M
= list(self
._multiplication
_table
) # copy
487 for i
in range(len(M
)):
488 # M had better be "square"
489 M
[i
] = [self
.monomial(i
)] + M
[i
]
490 M
= [["*"] + list(self
.gens())] + M
491 return table(M
, header_row
=True, header_column
=True, frame
=True)
494 def natural_basis(self
):
496 Return a more-natural representation of this algebra's basis.
498 Every finite-dimensional Euclidean Jordan Algebra is a direct
499 sum of five simple algebras, four of which comprise Hermitian
500 matrices. This method returns the original "natural" basis
501 for our underlying vector space. (Typically, the natural basis
502 is used to construct the multiplication table in the first place.)
504 Note that this will always return a matrix. The standard basis
505 in `R^n` will be returned as `n`-by-`1` column matrices.
509 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
510 ....: RealSymmetricEJA)
514 sage: J = RealSymmetricEJA(2)
516 Finite family {0: e0, 1: e1, 2: e2}
517 sage: J.natural_basis()
519 [1 0] [ 0 0.7071067811865475?] [0 0]
520 [0 0], [0.7071067811865475? 0], [0 1]
525 sage: J = JordanSpinEJA(2)
527 Finite family {0: e0, 1: e1}
528 sage: J.natural_basis()
535 if self
._natural
_basis
is None:
536 M
= self
.natural_basis_space()
537 return tuple( M(b
.to_vector()) for b
in self
.basis() )
539 return self
._natural
_basis
542 def natural_basis_space(self
):
544 Return the matrix space in which this algebra's natural basis
547 Generally this will be an `n`-by-`1` column-vector space,
548 except when the algebra is trivial. There it's `n`-by-`n`
549 (where `n` is zero), to ensure that two elements of the
550 natural basis space (empty matrices) can be multiplied.
552 if self
.is_trivial():
553 return MatrixSpace(self
.base_ring(), 0)
554 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
555 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
557 return self
._natural
_basis
[0].matrix_space()
563 Return the unit element of this algebra.
567 sage: from mjo.eja.eja_algebra import (HadamardEJA,
572 sage: J = HadamardEJA(5)
574 e0 + e1 + e2 + e3 + e4
578 The identity element acts like the identity::
580 sage: set_random_seed()
581 sage: J = random_eja()
582 sage: x = J.random_element()
583 sage: J.one()*x == x and x*J.one() == x
586 The matrix of the unit element's operator is the identity::
588 sage: set_random_seed()
589 sage: J = random_eja()
590 sage: actual = J.one().operator().matrix()
591 sage: expected = matrix.identity(J.base_ring(), J.dimension())
592 sage: actual == expected
595 Ensure that the cached unit element (often precomputed by
596 hand) agrees with the computed one::
598 sage: set_random_seed()
599 sage: J = random_eja()
600 sage: cached = J.one()
601 sage: J.one.clear_cache()
602 sage: J.one() == cached
606 # We can brute-force compute the matrices of the operators
607 # that correspond to the basis elements of this algebra.
608 # If some linear combination of those basis elements is the
609 # algebra identity, then the same linear combination of
610 # their matrices has to be the identity matrix.
612 # Of course, matrices aren't vectors in sage, so we have to
613 # appeal to the "long vectors" isometry.
614 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
616 # Now we use basic linear algebra to find the coefficients,
617 # of the matrices-as-vectors-linear-combination, which should
618 # work for the original algebra basis too.
619 A
= matrix(self
.base_ring(), oper_vecs
)
621 # We used the isometry on the left-hand side already, but we
622 # still need to do it for the right-hand side. Recall that we
623 # wanted something that summed to the identity matrix.
624 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
626 # Now if there's an identity element in the algebra, this
627 # should work. We solve on the left to avoid having to
628 # transpose the matrix "A".
629 return self
.from_vector(A
.solve_left(b
))
632 def peirce_decomposition(self
, c
):
634 The Peirce decomposition of this algebra relative to the
637 In the future, this can be extended to a complete system of
638 orthogonal idempotents.
642 - ``c`` -- an idempotent of this algebra.
646 A triple (J0, J5, J1) containing two subalgebras and one subspace
649 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
650 corresponding to the eigenvalue zero.
652 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
653 corresponding to the eigenvalue one-half.
655 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
656 corresponding to the eigenvalue one.
658 These are the only possible eigenspaces for that operator, and this
659 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
660 orthogonal, and are subalgebras of this algebra with the appropriate
665 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
669 The canonical example comes from the symmetric matrices, which
670 decompose into diagonal and off-diagonal parts::
672 sage: J = RealSymmetricEJA(3)
673 sage: C = matrix(QQ, [ [1,0,0],
677 sage: J0,J5,J1 = J.peirce_decomposition(c)
679 Euclidean Jordan algebra of dimension 1...
681 Vector space of degree 6 and dimension 2...
683 Euclidean Jordan algebra of dimension 3...
684 sage: J0.one().natural_representation()
688 sage: orig_df = AA.options.display_format
689 sage: AA.options.display_format = 'radical'
690 sage: J.from_vector(J5.basis()[0]).natural_representation()
694 sage: J.from_vector(J5.basis()[1]).natural_representation()
698 sage: AA.options.display_format = orig_df
699 sage: J1.one().natural_representation()
706 Every algebra decomposes trivially with respect to its identity
709 sage: set_random_seed()
710 sage: J = random_eja()
711 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
712 sage: J0.dimension() == 0 and J5.dimension() == 0
714 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
717 The decomposition is into eigenspaces, and its components are
718 therefore necessarily orthogonal. Moreover, the identity
719 elements in the two subalgebras are the projections onto their
720 respective subspaces of the superalgebra's identity element::
722 sage: set_random_seed()
723 sage: J = random_eja()
724 sage: x = J.random_element()
725 sage: if not J.is_trivial():
726 ....: while x.is_nilpotent():
727 ....: x = J.random_element()
728 sage: c = x.subalgebra_idempotent()
729 sage: J0,J5,J1 = J.peirce_decomposition(c)
731 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
732 ....: w = w.superalgebra_element()
733 ....: y = J.from_vector(y)
734 ....: z = z.superalgebra_element()
735 ....: ipsum += w.inner_product(y).abs()
736 ....: ipsum += w.inner_product(z).abs()
737 ....: ipsum += y.inner_product(z).abs()
740 sage: J1(c) == J1.one()
742 sage: J0(J.one() - c) == J0.one()
746 if not c
.is_idempotent():
747 raise ValueError("element is not idempotent: %s" % c
)
749 # Default these to what they should be if they turn out to be
750 # trivial, because eigenspaces_left() won't return eigenvalues
751 # corresponding to trivial spaces (e.g. it returns only the
752 # eigenspace corresponding to lambda=1 if you take the
753 # decomposition relative to the identity element).
754 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
755 J0
= trivial
# eigenvalue zero
756 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
757 J1
= trivial
# eigenvalue one
759 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
760 if eigval
== ~
(self
.base_ring()(2)):
763 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
764 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
772 raise ValueError("unexpected eigenvalue: %s" % eigval
)
777 def random_element(self
, thorough
=False):
779 Return a random element of this algebra.
781 Our algebra superclass method only returns a linear
782 combination of at most two basis elements. We instead
783 want the vector space "random element" method that
784 returns a more diverse selection.
788 - ``thorough`` -- (boolean; default False) whether or not we
789 should generate irrational coefficients for the random
790 element when our base ring is irrational; this slows the
791 algebra operations to a crawl, but any truly random method
795 # For a general base ring... maybe we can trust this to do the
796 # right thing? Unlikely, but.
797 V
= self
.vector_space()
798 v
= V
.random_element()
800 if self
.base_ring() is AA
:
801 # The "random element" method of the algebraic reals is
802 # stupid at the moment, and only returns integers between
803 # -2 and 2, inclusive:
805 # https://trac.sagemath.org/ticket/30875
807 # Instead, we implement our own "random vector" method,
808 # and then coerce that into the algebra. We use the vector
809 # space degree here instead of the dimension because a
810 # subalgebra could (for example) be spanned by only two
811 # vectors, each with five coordinates. We need to
812 # generate all five coordinates.
814 v
*= QQbar
.random_element().real()
816 v
*= QQ
.random_element()
818 return self
.from_vector(V
.coordinate_vector(v
))
820 def random_elements(self
, count
, thorough
=False):
822 Return ``count`` random elements as a tuple.
826 - ``thorough`` -- (boolean; default False) whether or not we
827 should generate irrational coefficients for the random
828 elements when our base ring is irrational; this slows the
829 algebra operations to a crawl, but any truly random method
834 sage: from mjo.eja.eja_algebra import JordanSpinEJA
838 sage: J = JordanSpinEJA(3)
839 sage: x,y,z = J.random_elements(3)
840 sage: all( [ x in J, y in J, z in J ])
842 sage: len( J.random_elements(10) ) == 10
846 return tuple( self
.random_element(thorough
)
847 for idx
in range(count
) )
851 def _charpoly_coefficients(self
):
853 The `r` polynomial coefficients of the "characteristic polynomial
857 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
858 R
= PolynomialRing(self
.base_ring(), var_names
)
860 F
= R
.fraction_field()
863 # From a result in my book, these are the entries of the
864 # basis representation of L_x.
865 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
868 L_x
= matrix(F
, n
, n
, L_x_i_j
)
871 if self
.rank
.is_in_cache():
873 # There's no need to pad the system with redundant
874 # columns if we *know* they'll be redundant.
877 # Compute an extra power in case the rank is equal to
878 # the dimension (otherwise, we would stop at x^(r-1)).
879 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
880 for k
in range(n
+1) ]
881 A
= matrix
.column(F
, x_powers
[:n
])
882 AE
= A
.extended_echelon_form()
889 # The theory says that only the first "r" coefficients are
890 # nonzero, and they actually live in the original polynomial
891 # ring and not the fraction field. We negate them because
892 # in the actual characteristic polynomial, they get moved
893 # to the other side where x^r lives.
894 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
899 Return the rank of this EJA.
901 This is a cached method because we know the rank a priori for
902 all of the algebras we can construct. Thus we can avoid the
903 expensive ``_charpoly_coefficients()`` call unless we truly
904 need to compute the whole characteristic polynomial.
908 sage: from mjo.eja.eja_algebra import (HadamardEJA,
910 ....: RealSymmetricEJA,
911 ....: ComplexHermitianEJA,
912 ....: QuaternionHermitianEJA,
917 The rank of the Jordan spin algebra is always two::
919 sage: JordanSpinEJA(2).rank()
921 sage: JordanSpinEJA(3).rank()
923 sage: JordanSpinEJA(4).rank()
926 The rank of the `n`-by-`n` Hermitian real, complex, or
927 quaternion matrices is `n`::
929 sage: RealSymmetricEJA(4).rank()
931 sage: ComplexHermitianEJA(3).rank()
933 sage: QuaternionHermitianEJA(2).rank()
938 Ensure that every EJA that we know how to construct has a
939 positive integer rank, unless the algebra is trivial in
940 which case its rank will be zero::
942 sage: set_random_seed()
943 sage: J = random_eja()
947 sage: r > 0 or (r == 0 and J.is_trivial())
950 Ensure that computing the rank actually works, since the ranks
951 of all simple algebras are known and will be cached by default::
953 sage: J = HadamardEJA(4)
954 sage: J.rank.clear_cache()
960 sage: J = JordanSpinEJA(4)
961 sage: J.rank.clear_cache()
967 sage: J = RealSymmetricEJA(3)
968 sage: J.rank.clear_cache()
974 sage: J = ComplexHermitianEJA(2)
975 sage: J.rank.clear_cache()
981 sage: J = QuaternionHermitianEJA(2)
982 sage: J.rank.clear_cache()
986 return len(self
._charpoly
_coefficients
())
989 def vector_space(self
):
991 Return the vector space that underlies this algebra.
995 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
999 sage: J = RealSymmetricEJA(2)
1000 sage: J.vector_space()
1001 Vector space of dimension 3 over...
1004 return self
.zero().to_vector().parent().ambient_vector_space()
1007 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1010 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1012 Algebras whose basis consists of vectors with rational
1013 entries. Equivalently, algebras whose multiplication tables
1014 contain only rational coefficients.
1016 When an EJA has a basis that can be made rational, we can speed up
1017 the computation of its characteristic polynomial by doing it over
1018 ``QQ``. All of the named EJA constructors that we provide fall
1022 def _charpoly_coefficients(self
):
1024 Override the parent method with something that tries to compute
1025 over a faster (non-extension) field.
1029 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1033 The base ring of the resulting polynomial coefficients is what
1034 it should be, and not the rationals (unless the algebra was
1035 already over the rationals)::
1037 sage: J = JordanSpinEJA(3)
1038 sage: J._charpoly_coefficients()
1039 (X1^2 - X2^2 - X3^2, -2*X1)
1040 sage: a0 = J._charpoly_coefficients()[0]
1042 Algebraic Real Field
1043 sage: a0.base_ring()
1044 Algebraic Real Field
1047 if self
.base_ring() is QQ
:
1048 # There's no need to construct *another* algebra over the
1049 # rationals if this one is already over the rationals.
1050 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1051 return superclass
._charpoly
_coefficients
()
1054 map(lambda x
: x
.to_vector(), ls
)
1055 for ls
in self
._multiplication
_table
1058 # Do the computation over the rationals. The answer will be
1059 # the same, because our basis coordinates are (essentially)
1061 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1065 a
= J
._charpoly
_coefficients
()
1066 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1069 class ConcreteEuclideanJordanAlgebra
:
1071 A class for the Euclidean Jordan algebras that we know by name.
1073 These are the Jordan algebras whose basis, multiplication table,
1074 rank, and so on are known a priori. More to the point, they are
1075 the Euclidean Jordan algebras for which we are able to conjure up
1076 a "random instance."
1080 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1084 Our natural basis is normalized with respect to the natural inner
1085 product unless we specify otherwise::
1087 sage: set_random_seed()
1088 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1089 sage: all( b.norm() == 1 for b in J.gens() )
1092 Since our natural basis is normalized with respect to the natural
1093 inner product, and since we know that this algebra is an EJA, any
1094 left-multiplication operator's matrix will be symmetric because
1095 natural->EJA basis representation is an isometry and within the EJA
1096 the operator is self-adjoint by the Jordan axiom::
1098 sage: set_random_seed()
1099 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1100 sage: x = J.random_element()
1101 sage: x.operator().matrix().is_symmetric()
1107 def _max_random_instance_size():
1109 Return an integer "size" that is an upper bound on the size of
1110 this algebra when it is used in a random test
1111 case. Unfortunately, the term "size" is ambiguous -- when
1112 dealing with `R^n` under either the Hadamard or Jordan spin
1113 product, the "size" refers to the dimension `n`. When dealing
1114 with a matrix algebra (real symmetric or complex/quaternion
1115 Hermitian), it refers to the size of the matrix, which is far
1116 less than the dimension of the underlying vector space.
1118 This method must be implemented in each subclass.
1120 raise NotImplementedError
1123 def random_instance(cls
, field
=AA
, **kwargs
):
1125 Return a random instance of this type of algebra.
1127 This method should be implemented in each subclass.
1129 from sage
.misc
.prandom
import choice
1130 eja_class
= choice(cls
.__subclasses
__())
1131 return eja_class
.random_instance(field
)
1134 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1136 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1138 Compared to the superclass constructor, we take a basis instead of
1139 a multiplication table because the latter can be computed in terms
1140 of the former when the product is known (like it is here).
1142 # Used in this class's fast _charpoly_coefficients() override.
1143 self
._basis
_normalizers
= None
1145 # We're going to loop through this a few times, so now's a good
1146 # time to ensure that it isn't a generator expression.
1147 basis
= tuple(basis
)
1149 algebra_dim
= len(basis
)
1150 degree
= 0 # size of the matrices
1152 degree
= basis
[0].nrows()
1154 if algebra_dim
> 1 and normalize_basis
:
1155 # We'll need sqrt(2) to normalize the basis, and this
1156 # winds up in the multiplication table, so the whole
1157 # algebra needs to be over the field extension.
1158 R
= PolynomialRing(field
, 'z')
1161 if p
.is_irreducible():
1162 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1163 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1164 self
._basis
_normalizers
= tuple(
1165 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1166 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1168 # Now compute the multiplication and inner product tables.
1169 # We have to do this *after* normalizing the basis, because
1170 # scaling affects the answers.
1171 V
= VectorSpace(field
, degree
**2)
1172 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1173 mult_table
= [[W
.zero() for j
in range(algebra_dim
)]
1174 for i
in range(algebra_dim
)]
1175 ip_table
= [[W
.zero() for j
in range(algebra_dim
)]
1176 for i
in range(algebra_dim
)]
1177 for i
in range(algebra_dim
):
1178 for j
in range(algebra_dim
):
1179 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1180 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1183 # HACK: ignore the error here if we don't need the
1184 # inner product (as is the case when we construct
1185 # a dummy QQ-algebra for fast charpoly coefficients.
1186 ip_table
[i
][j
] = self
.natural_inner_product(basis
[i
],
1193 self
._inner
_product
_matrix
= matrix(field
,ip_table
)
1197 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1199 natural_basis
=basis
,
1202 if algebra_dim
== 0:
1203 self
.one
.set_cache(self
.zero())
1205 n
= basis
[0].nrows()
1206 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1207 # details of this algebra.
1208 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1212 def _charpoly_coefficients(self
):
1214 Override the parent method with something that tries to compute
1215 over a faster (non-extension) field.
1217 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1218 # We didn't normalize, or the basis we started with had
1219 # entries in a nice field already. Just compute the thing.
1220 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1222 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1223 self
._basis
_normalizers
) )
1225 # Do this over the rationals and convert back at the end.
1226 # Only works because we know the entries of the basis are
1227 # integers. The argument ``check_axioms=False`` is required
1228 # because the trace inner-product method for this
1229 # class is a stub and can't actually be checked.
1230 J
= MatrixEuclideanJordanAlgebra(QQ
,
1232 normalize_basis
=False,
1235 a
= J
._charpoly
_coefficients
()
1237 # Unfortunately, changing the basis does change the
1238 # coefficients of the characteristic polynomial, but since
1239 # these are really the coefficients of the "characteristic
1240 # polynomial of" function, everything is still nice and
1241 # unevaluated. It's therefore "obvious" how scaling the
1242 # basis affects the coordinate variables X1, X2, et
1243 # cetera. Scaling the first basis vector up by "n" adds a
1244 # factor of 1/n into every "X1" term, for example. So here
1245 # we simply undo the basis_normalizer scaling that we
1246 # performed earlier.
1248 # The a[0] access here is safe because trivial algebras
1249 # won't have any basis normalizers and therefore won't
1250 # make it to this "else" branch.
1251 XS
= a
[0].parent().gens()
1252 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1253 for i
in range(len(XS
)) }
1254 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1260 Embed the matrix ``M`` into a space of real matrices.
1262 The matrix ``M`` can have entries in any field at the moment:
1263 the real numbers, complex numbers, or quaternions. And although
1264 they are not a field, we can probably support octonions at some
1265 point, too. This function returns a real matrix that "acts like"
1266 the original with respect to matrix multiplication; i.e.
1268 real_embed(M*N) = real_embed(M)*real_embed(N)
1271 raise NotImplementedError
1275 def real_unembed(M
):
1277 The inverse of :meth:`real_embed`.
1279 raise NotImplementedError
1282 def natural_inner_product(cls
,X
,Y
):
1283 Xu
= cls
.real_unembed(X
)
1284 Yu
= cls
.real_unembed(Y
)
1285 tr
= (Xu
*Yu
).trace()
1288 # Works in QQ, AA, RDF, et cetera.
1290 except AttributeError:
1291 # A quaternion doesn't have a real() method, but does
1292 # have coefficient_tuple() method that returns the
1293 # coefficients of 1, i, j, and k -- in that order.
1294 return tr
.coefficient_tuple()[0]
1297 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1301 The identity function, for embedding real matrices into real
1307 def real_unembed(M
):
1309 The identity function, for unembedding real matrices from real
1315 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1316 ConcreteEuclideanJordanAlgebra
):
1318 The rank-n simple EJA consisting of real symmetric n-by-n
1319 matrices, the usual symmetric Jordan product, and the trace inner
1320 product. It has dimension `(n^2 + n)/2` over the reals.
1324 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1328 sage: J = RealSymmetricEJA(2)
1329 sage: e0, e1, e2 = J.gens()
1337 In theory, our "field" can be any subfield of the reals::
1339 sage: RealSymmetricEJA(2, RDF)
1340 Euclidean Jordan algebra of dimension 3 over Real Double Field
1341 sage: RealSymmetricEJA(2, RR)
1342 Euclidean Jordan algebra of dimension 3 over Real Field with
1343 53 bits of precision
1347 The dimension of this algebra is `(n^2 + n) / 2`::
1349 sage: set_random_seed()
1350 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1351 sage: n = ZZ.random_element(1, n_max)
1352 sage: J = RealSymmetricEJA(n)
1353 sage: J.dimension() == (n^2 + n)/2
1356 The Jordan multiplication is what we think it is::
1358 sage: set_random_seed()
1359 sage: J = RealSymmetricEJA.random_instance()
1360 sage: x,y = J.random_elements(2)
1361 sage: actual = (x*y).natural_representation()
1362 sage: X = x.natural_representation()
1363 sage: Y = y.natural_representation()
1364 sage: expected = (X*Y + Y*X)/2
1365 sage: actual == expected
1367 sage: J(expected) == x*y
1370 We can change the generator prefix::
1372 sage: RealSymmetricEJA(3, prefix='q').gens()
1373 (q0, q1, q2, q3, q4, q5)
1375 We can construct the (trivial) algebra of rank zero::
1377 sage: RealSymmetricEJA(0)
1378 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1382 def _denormalized_basis(cls
, n
, field
):
1384 Return a basis for the space of real symmetric n-by-n matrices.
1388 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1392 sage: set_random_seed()
1393 sage: n = ZZ.random_element(1,5)
1394 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1395 sage: all( M.is_symmetric() for M in B)
1399 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1403 for j
in range(i
+1):
1404 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1408 Sij
= Eij
+ Eij
.transpose()
1414 def _max_random_instance_size():
1415 return 4 # Dimension 10
1418 def random_instance(cls
, field
=AA
, **kwargs
):
1420 Return a random instance of this type of algebra.
1422 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1423 return cls(n
, field
, **kwargs
)
1425 def __init__(self
, n
, field
=AA
, **kwargs
):
1426 basis
= self
._denormalized
_basis
(n
, field
)
1427 super(RealSymmetricEJA
, self
).__init
__(field
,
1431 self
.rank
.set_cache(n
)
1434 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1438 Embed the n-by-n complex matrix ``M`` into the space of real
1439 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1440 bi` to the block matrix ``[[a,b],[-b,a]]``.
1444 sage: from mjo.eja.eja_algebra import \
1445 ....: ComplexMatrixEuclideanJordanAlgebra
1449 sage: F = QuadraticField(-1, 'I')
1450 sage: x1 = F(4 - 2*i)
1451 sage: x2 = F(1 + 2*i)
1454 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1455 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1464 Embedding is a homomorphism (isomorphism, in fact)::
1466 sage: set_random_seed()
1467 sage: n = ZZ.random_element(3)
1468 sage: F = QuadraticField(-1, 'I')
1469 sage: X = random_matrix(F, n)
1470 sage: Y = random_matrix(F, n)
1471 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1472 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1473 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1480 raise ValueError("the matrix 'M' must be square")
1482 # We don't need any adjoined elements...
1483 field
= M
.base_ring().base_ring()
1487 a
= z
.list()[0] # real part, I guess
1488 b
= z
.list()[1] # imag part, I guess
1489 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1491 return matrix
.block(field
, n
, blocks
)
1495 def real_unembed(M
):
1497 The inverse of _embed_complex_matrix().
1501 sage: from mjo.eja.eja_algebra import \
1502 ....: ComplexMatrixEuclideanJordanAlgebra
1506 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1507 ....: [-2, 1, -4, 3],
1508 ....: [ 9, 10, 11, 12],
1509 ....: [-10, 9, -12, 11] ])
1510 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1512 [ 10*I + 9 12*I + 11]
1516 Unembedding is the inverse of embedding::
1518 sage: set_random_seed()
1519 sage: F = QuadraticField(-1, 'I')
1520 sage: M = random_matrix(F, 3)
1521 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1522 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1528 raise ValueError("the matrix 'M' must be square")
1529 if not n
.mod(2).is_zero():
1530 raise ValueError("the matrix 'M' must be a complex embedding")
1532 # If "M" was normalized, its base ring might have roots
1533 # adjoined and they can stick around after unembedding.
1534 field
= M
.base_ring()
1535 R
= PolynomialRing(field
, 'z')
1538 # Sage doesn't know how to embed AA into QQbar, i.e. how
1539 # to adjoin sqrt(-1) to AA.
1542 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1545 # Go top-left to bottom-right (reading order), converting every
1546 # 2-by-2 block we see to a single complex element.
1548 for k
in range(n
/2):
1549 for j
in range(n
/2):
1550 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1551 if submat
[0,0] != submat
[1,1]:
1552 raise ValueError('bad on-diagonal submatrix')
1553 if submat
[0,1] != -submat
[1,0]:
1554 raise ValueError('bad off-diagonal submatrix')
1555 z
= submat
[0,0] + submat
[0,1]*i
1558 return matrix(F
, n
/2, elements
)
1562 def natural_inner_product(cls
,X
,Y
):
1564 Compute a natural inner product in this algebra directly from
1569 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1573 This gives the same answer as the slow, default method implemented
1574 in :class:`MatrixEuclideanJordanAlgebra`::
1576 sage: set_random_seed()
1577 sage: J = ComplexHermitianEJA.random_instance()
1578 sage: x,y = J.random_elements(2)
1579 sage: Xe = x.natural_representation()
1580 sage: Ye = y.natural_representation()
1581 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1582 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1583 sage: expected = (X*Y).trace().real()
1584 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1585 sage: actual == expected
1589 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1592 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1593 ConcreteEuclideanJordanAlgebra
):
1595 The rank-n simple EJA consisting of complex Hermitian n-by-n
1596 matrices over the real numbers, the usual symmetric Jordan product,
1597 and the real-part-of-trace inner product. It has dimension `n^2` over
1602 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1606 In theory, our "field" can be any subfield of the reals::
1608 sage: ComplexHermitianEJA(2, RDF)
1609 Euclidean Jordan algebra of dimension 4 over Real Double Field
1610 sage: ComplexHermitianEJA(2, RR)
1611 Euclidean Jordan algebra of dimension 4 over Real Field with
1612 53 bits of precision
1616 The dimension of this algebra is `n^2`::
1618 sage: set_random_seed()
1619 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1620 sage: n = ZZ.random_element(1, n_max)
1621 sage: J = ComplexHermitianEJA(n)
1622 sage: J.dimension() == n^2
1625 The Jordan multiplication is what we think it is::
1627 sage: set_random_seed()
1628 sage: J = ComplexHermitianEJA.random_instance()
1629 sage: x,y = J.random_elements(2)
1630 sage: actual = (x*y).natural_representation()
1631 sage: X = x.natural_representation()
1632 sage: Y = y.natural_representation()
1633 sage: expected = (X*Y + Y*X)/2
1634 sage: actual == expected
1636 sage: J(expected) == x*y
1639 We can change the generator prefix::
1641 sage: ComplexHermitianEJA(2, prefix='z').gens()
1644 We can construct the (trivial) algebra of rank zero::
1646 sage: ComplexHermitianEJA(0)
1647 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1652 def _denormalized_basis(cls
, n
, field
):
1654 Returns a basis for the space of complex Hermitian n-by-n matrices.
1656 Why do we embed these? Basically, because all of numerical linear
1657 algebra assumes that you're working with vectors consisting of `n`
1658 entries from a field and scalars from the same field. There's no way
1659 to tell SageMath that (for example) the vectors contain complex
1660 numbers, while the scalar field is real.
1664 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1668 sage: set_random_seed()
1669 sage: n = ZZ.random_element(1,5)
1670 sage: field = QuadraticField(2, 'sqrt2')
1671 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1672 sage: all( M.is_symmetric() for M in B)
1676 R
= PolynomialRing(field
, 'z')
1678 F
= field
.extension(z
**2 + 1, 'I')
1681 # This is like the symmetric case, but we need to be careful:
1683 # * We want conjugate-symmetry, not just symmetry.
1684 # * The diagonal will (as a result) be real.
1688 for j
in range(i
+1):
1689 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1691 Sij
= cls
.real_embed(Eij
)
1694 # The second one has a minus because it's conjugated.
1695 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1697 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1700 # Since we embedded these, we can drop back to the "field" that we
1701 # started with instead of the complex extension "F".
1702 return ( s
.change_ring(field
) for s
in S
)
1705 def __init__(self
, n
, field
=AA
, **kwargs
):
1706 basis
= self
._denormalized
_basis
(n
,field
)
1707 super(ComplexHermitianEJA
,self
).__init
__(field
,
1711 self
.rank
.set_cache(n
)
1714 def _max_random_instance_size():
1715 return 3 # Dimension 9
1718 def random_instance(cls
, field
=AA
, **kwargs
):
1720 Return a random instance of this type of algebra.
1722 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1723 return cls(n
, field
, **kwargs
)
1725 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1729 Embed the n-by-n quaternion matrix ``M`` into the space of real
1730 matrices of size 4n-by-4n by first sending each quaternion entry `z
1731 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1732 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1737 sage: from mjo.eja.eja_algebra import \
1738 ....: QuaternionMatrixEuclideanJordanAlgebra
1742 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1743 sage: i,j,k = Q.gens()
1744 sage: x = 1 + 2*i + 3*j + 4*k
1745 sage: M = matrix(Q, 1, [[x]])
1746 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1752 Embedding is a homomorphism (isomorphism, in fact)::
1754 sage: set_random_seed()
1755 sage: n = ZZ.random_element(2)
1756 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1757 sage: X = random_matrix(Q, n)
1758 sage: Y = random_matrix(Q, n)
1759 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1760 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1761 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1766 quaternions
= M
.base_ring()
1769 raise ValueError("the matrix 'M' must be square")
1771 F
= QuadraticField(-1, 'I')
1776 t
= z
.coefficient_tuple()
1781 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1782 [-c
+ d
*i
, a
- b
*i
]])
1783 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1784 blocks
.append(realM
)
1786 # We should have real entries by now, so use the realest field
1787 # we've got for the return value.
1788 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1793 def real_unembed(M
):
1795 The inverse of _embed_quaternion_matrix().
1799 sage: from mjo.eja.eja_algebra import \
1800 ....: QuaternionMatrixEuclideanJordanAlgebra
1804 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1805 ....: [-2, 1, -4, 3],
1806 ....: [-3, 4, 1, -2],
1807 ....: [-4, -3, 2, 1]])
1808 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1809 [1 + 2*i + 3*j + 4*k]
1813 Unembedding is the inverse of embedding::
1815 sage: set_random_seed()
1816 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1817 sage: M = random_matrix(Q, 3)
1818 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1819 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1825 raise ValueError("the matrix 'M' must be square")
1826 if not n
.mod(4).is_zero():
1827 raise ValueError("the matrix 'M' must be a quaternion embedding")
1829 # Use the base ring of the matrix to ensure that its entries can be
1830 # multiplied by elements of the quaternion algebra.
1831 field
= M
.base_ring()
1832 Q
= QuaternionAlgebra(field
,-1,-1)
1835 # Go top-left to bottom-right (reading order), converting every
1836 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1839 for l
in range(n
/4):
1840 for m
in range(n
/4):
1841 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1842 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1843 if submat
[0,0] != submat
[1,1].conjugate():
1844 raise ValueError('bad on-diagonal submatrix')
1845 if submat
[0,1] != -submat
[1,0].conjugate():
1846 raise ValueError('bad off-diagonal submatrix')
1847 z
= submat
[0,0].real()
1848 z
+= submat
[0,0].imag()*i
1849 z
+= submat
[0,1].real()*j
1850 z
+= submat
[0,1].imag()*k
1853 return matrix(Q
, n
/4, elements
)
1857 def natural_inner_product(cls
,X
,Y
):
1859 Compute a natural inner product in this algebra directly from
1864 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1868 This gives the same answer as the slow, default method implemented
1869 in :class:`MatrixEuclideanJordanAlgebra`::
1871 sage: set_random_seed()
1872 sage: J = QuaternionHermitianEJA.random_instance()
1873 sage: x,y = J.random_elements(2)
1874 sage: Xe = x.natural_representation()
1875 sage: Ye = y.natural_representation()
1876 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1877 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1878 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1879 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1880 sage: actual == expected
1884 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1887 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1888 ConcreteEuclideanJordanAlgebra
):
1890 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1891 matrices, the usual symmetric Jordan product, and the
1892 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1897 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1901 In theory, our "field" can be any subfield of the reals::
1903 sage: QuaternionHermitianEJA(2, RDF)
1904 Euclidean Jordan algebra of dimension 6 over Real Double Field
1905 sage: QuaternionHermitianEJA(2, RR)
1906 Euclidean Jordan algebra of dimension 6 over Real Field with
1907 53 bits of precision
1911 The dimension of this algebra is `2*n^2 - n`::
1913 sage: set_random_seed()
1914 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
1915 sage: n = ZZ.random_element(1, n_max)
1916 sage: J = QuaternionHermitianEJA(n)
1917 sage: J.dimension() == 2*(n^2) - n
1920 The Jordan multiplication is what we think it is::
1922 sage: set_random_seed()
1923 sage: J = QuaternionHermitianEJA.random_instance()
1924 sage: x,y = J.random_elements(2)
1925 sage: actual = (x*y).natural_representation()
1926 sage: X = x.natural_representation()
1927 sage: Y = y.natural_representation()
1928 sage: expected = (X*Y + Y*X)/2
1929 sage: actual == expected
1931 sage: J(expected) == x*y
1934 We can change the generator prefix::
1936 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1937 (a0, a1, a2, a3, a4, a5)
1939 We can construct the (trivial) algebra of rank zero::
1941 sage: QuaternionHermitianEJA(0)
1942 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1946 def _denormalized_basis(cls
, n
, field
):
1948 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1950 Why do we embed these? Basically, because all of numerical
1951 linear algebra assumes that you're working with vectors consisting
1952 of `n` entries from a field and scalars from the same field. There's
1953 no way to tell SageMath that (for example) the vectors contain
1954 complex numbers, while the scalar field is real.
1958 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1962 sage: set_random_seed()
1963 sage: n = ZZ.random_element(1,5)
1964 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1965 sage: all( M.is_symmetric() for M in B )
1969 Q
= QuaternionAlgebra(QQ
,-1,-1)
1972 # This is like the symmetric case, but we need to be careful:
1974 # * We want conjugate-symmetry, not just symmetry.
1975 # * The diagonal will (as a result) be real.
1979 for j
in range(i
+1):
1980 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1982 Sij
= cls
.real_embed(Eij
)
1985 # The second, third, and fourth ones have a minus
1986 # because they're conjugated.
1987 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1989 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1991 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1993 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1996 # Since we embedded these, we can drop back to the "field" that we
1997 # started with instead of the quaternion algebra "Q".
1998 return ( s
.change_ring(field
) for s
in S
)
2001 def __init__(self
, n
, field
=AA
, **kwargs
):
2002 basis
= self
._denormalized
_basis
(n
,field
)
2003 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2007 self
.rank
.set_cache(n
)
2010 def _max_random_instance_size():
2012 The maximum rank of a random QuaternionHermitianEJA.
2014 return 2 # Dimension 6
2017 def random_instance(cls
, field
=AA
, **kwargs
):
2019 Return a random instance of this type of algebra.
2021 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2022 return cls(n
, field
, **kwargs
)
2025 class HadamardEJA(RationalBasisEuclideanJordanAlgebra
,
2026 ConcreteEuclideanJordanAlgebra
):
2028 Return the Euclidean Jordan Algebra corresponding to the set
2029 `R^n` under the Hadamard product.
2031 Note: this is nothing more than the Cartesian product of ``n``
2032 copies of the spin algebra. Once Cartesian product algebras
2033 are implemented, this can go.
2037 sage: from mjo.eja.eja_algebra import HadamardEJA
2041 This multiplication table can be verified by hand::
2043 sage: J = HadamardEJA(3)
2044 sage: e0,e1,e2 = J.gens()
2060 We can change the generator prefix::
2062 sage: HadamardEJA(3, prefix='r').gens()
2066 def __init__(self
, n
, field
=AA
, **kwargs
):
2067 V
= VectorSpace(field
, n
)
2068 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
2071 # Inner products are real numbers and not algebra
2072 # elements, so once we turn the algebra element
2073 # into a vector in inner_product(), we never go
2074 # back. As a result -- contrary to what we do with
2075 # self._multiplication_table -- we store the inner
2076 # product table as a plain old matrix and not as
2077 # an algebra operator.
2078 ip_table
= matrix
.identity(field
,n
)
2079 self
._inner
_product
_matrix
= ip_table
2081 super(HadamardEJA
, self
).__init
__(field
,
2085 self
.rank
.set_cache(n
)
2088 self
.one
.set_cache( self
.zero() )
2090 self
.one
.set_cache( sum(self
.gens()) )
2093 def _max_random_instance_size():
2095 The maximum dimension of a random HadamardEJA.
2100 def random_instance(cls
, field
=AA
, **kwargs
):
2102 Return a random instance of this type of algebra.
2104 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2105 return cls(n
, field
, **kwargs
)
2108 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
,
2109 ConcreteEuclideanJordanAlgebra
):
2111 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2112 with the half-trace inner product and jordan product ``x*y =
2113 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2114 a symmetric positive-definite "bilinear form" matrix. Its
2115 dimension is the size of `B`, and it has rank two in dimensions
2116 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2117 the identity matrix of order ``n``.
2119 We insist that the one-by-one upper-left identity block of `B` be
2120 passed in as well so that we can be passed a matrix of size zero
2121 to construct a trivial algebra.
2125 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2126 ....: JordanSpinEJA)
2130 When no bilinear form is specified, the identity matrix is used,
2131 and the resulting algebra is the Jordan spin algebra::
2133 sage: B = matrix.identity(AA,3)
2134 sage: J0 = BilinearFormEJA(B)
2135 sage: J1 = JordanSpinEJA(3)
2136 sage: J0.multiplication_table() == J0.multiplication_table()
2139 An error is raised if the matrix `B` does not correspond to a
2140 positive-definite bilinear form::
2142 sage: B = matrix.random(QQ,2,3)
2143 sage: J = BilinearFormEJA(B)
2144 Traceback (most recent call last):
2146 ValueError: bilinear form is not positive-definite
2147 sage: B = matrix.zero(QQ,3)
2148 sage: J = BilinearFormEJA(B)
2149 Traceback (most recent call last):
2151 ValueError: bilinear form is not positive-definite
2155 We can create a zero-dimensional algebra::
2157 sage: B = matrix.identity(AA,0)
2158 sage: J = BilinearFormEJA(B)
2162 We can check the multiplication condition given in the Jordan, von
2163 Neumann, and Wigner paper (and also discussed on my "On the
2164 symmetry..." paper). Note that this relies heavily on the standard
2165 choice of basis, as does anything utilizing the bilinear form matrix::
2167 sage: set_random_seed()
2168 sage: n = ZZ.random_element(5)
2169 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2170 sage: B11 = matrix.identity(QQ,1)
2171 sage: B22 = M.transpose()*M
2172 sage: B = block_matrix(2,2,[ [B11,0 ],
2174 sage: J = BilinearFormEJA(B)
2175 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2176 sage: V = J.vector_space()
2177 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2178 ....: for ei in eis ]
2179 sage: actual = [ sis[i]*sis[j]
2180 ....: for i in range(n-1)
2181 ....: for j in range(n-1) ]
2182 sage: expected = [ J.one() if i == j else J.zero()
2183 ....: for i in range(n-1)
2184 ....: for j in range(n-1) ]
2185 sage: actual == expected
2188 def __init__(self
, B
, field
=AA
, **kwargs
):
2191 if not B
.is_positive_definite():
2192 raise ValueError("bilinear form is not positive-definite")
2194 V
= VectorSpace(field
, n
)
2195 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2204 z0
= (B
*x
).inner_product(y
)
2205 zbar
= y0
*xbar
+ x0
*ybar
2206 z
= V([z0
] + zbar
.list())
2207 mult_table
[i
][j
] = z
2209 # Inner products are real numbers and not algebra
2210 # elements, so once we turn the algebra element
2211 # into a vector in inner_product(), we never go
2212 # back. As a result -- contrary to what we do with
2213 # self._multiplication_table -- we store the inner
2214 # product table as a plain old matrix and not as
2215 # an algebra operator.
2217 self
._inner
_product
_matrix
= ip_table
2219 super(BilinearFormEJA
, self
).__init
__(field
,
2224 # The rank of this algebra is two, unless we're in a
2225 # one-dimensional ambient space (because the rank is bounded
2226 # by the ambient dimension).
2227 self
.rank
.set_cache(min(n
,2))
2230 self
.one
.set_cache( self
.zero() )
2232 self
.one
.set_cache( self
.monomial(0) )
2235 def _max_random_instance_size():
2237 The maximum dimension of a random BilinearFormEJA.
2242 def random_instance(cls
, field
=AA
, **kwargs
):
2244 Return a random instance of this algebra.
2246 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2248 B
= matrix
.identity(field
, n
)
2249 return cls(B
, field
, **kwargs
)
2251 B11
= matrix
.identity(field
,1)
2252 M
= matrix
.random(field
, n
-1)
2253 I
= matrix
.identity(field
, n
-1)
2254 alpha
= field
.zero()
2255 while alpha
.is_zero():
2256 alpha
= field
.random_element().abs()
2257 B22
= M
.transpose()*M
+ alpha
*I
2259 from sage
.matrix
.special
import block_matrix
2260 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2263 return cls(B
, field
, **kwargs
)
2266 class JordanSpinEJA(BilinearFormEJA
):
2268 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2269 with the usual inner product and jordan product ``x*y =
2270 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2275 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2279 This multiplication table can be verified by hand::
2281 sage: J = JordanSpinEJA(4)
2282 sage: e0,e1,e2,e3 = J.gens()
2298 We can change the generator prefix::
2300 sage: JordanSpinEJA(2, prefix='B').gens()
2305 Ensure that we have the usual inner product on `R^n`::
2307 sage: set_random_seed()
2308 sage: J = JordanSpinEJA.random_instance()
2309 sage: x,y = J.random_elements(2)
2310 sage: actual = x.inner_product(y)
2311 sage: expected = x.to_vector().inner_product(y.to_vector())
2312 sage: actual == expected
2316 def __init__(self
, n
, field
=AA
, **kwargs
):
2317 # This is a special case of the BilinearFormEJA with the identity
2318 # matrix as its bilinear form.
2319 B
= matrix
.identity(field
, n
)
2320 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2323 def _max_random_instance_size():
2325 The maximum dimension of a random JordanSpinEJA.
2330 def random_instance(cls
, field
=AA
, **kwargs
):
2332 Return a random instance of this type of algebra.
2334 Needed here to override the implementation for ``BilinearFormEJA``.
2336 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2337 return cls(n
, field
, **kwargs
)
2340 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
,
2341 ConcreteEuclideanJordanAlgebra
):
2343 The trivial Euclidean Jordan algebra consisting of only a zero element.
2347 sage: from mjo.eja.eja_algebra import TrivialEJA
2351 sage: J = TrivialEJA()
2358 sage: 7*J.one()*12*J.one()
2360 sage: J.one().inner_product(J.one())
2362 sage: J.one().norm()
2364 sage: J.one().subalgebra_generated_by()
2365 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2370 def __init__(self
, field
=AA
, **kwargs
):
2372 self
._inner
_product
_matrix
= matrix(field
,0)
2373 super(TrivialEJA
, self
).__init
__(field
,
2377 # The rank is zero using my definition, namely the dimension of the
2378 # largest subalgebra generated by any element.
2379 self
.rank
.set_cache(0)
2380 self
.one
.set_cache( self
.zero() )
2383 def random_instance(cls
, field
=AA
, **kwargs
):
2384 # We don't take a "size" argument so the superclass method is
2385 # inappropriate for us.
2386 return cls(field
, **kwargs
)
2388 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2390 The external (orthogonal) direct sum of two other Euclidean Jordan
2391 algebras. Essentially the Cartesian product of its two factors.
2392 Every Euclidean Jordan algebra decomposes into an orthogonal
2393 direct sum of simple Euclidean Jordan algebras, so no generality
2394 is lost by providing only this construction.
2398 sage: from mjo.eja.eja_algebra import (random_eja,
2400 ....: RealSymmetricEJA,
2405 sage: J1 = HadamardEJA(2)
2406 sage: J2 = RealSymmetricEJA(3)
2407 sage: J = DirectSumEJA(J1,J2)
2415 The external direct sum construction is only valid when the two factors
2416 have the same base ring; an error is raised otherwise::
2418 sage: set_random_seed()
2419 sage: J1 = random_eja(AA)
2420 sage: J2 = random_eja(QQ)
2421 sage: J = DirectSumEJA(J1,J2)
2422 Traceback (most recent call last):
2424 ValueError: algebras must share the same base field
2427 def __init__(self
, J1
, J2
, **kwargs
):
2428 if J1
.base_ring() != J2
.base_ring():
2429 raise ValueError("algebras must share the same base field")
2430 field
= J1
.base_ring()
2432 self
._factors
= (J1
, J2
)
2436 V
= VectorSpace(field
, n
)
2437 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2441 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2442 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2446 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2447 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2449 super(DirectSumEJA
, self
).__init
__(field
,
2453 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2458 Return the pair of this algebra's factors.
2462 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2463 ....: JordanSpinEJA,
2468 sage: J1 = HadamardEJA(2,QQ)
2469 sage: J2 = JordanSpinEJA(3,QQ)
2470 sage: J = DirectSumEJA(J1,J2)
2472 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2473 Euclidean Jordan algebra of dimension 3 over Rational Field)
2476 return self
._factors
2478 def projections(self
):
2480 Return a pair of projections onto this algebra's factors.
2484 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2485 ....: ComplexHermitianEJA,
2490 sage: J1 = JordanSpinEJA(2)
2491 sage: J2 = ComplexHermitianEJA(2)
2492 sage: J = DirectSumEJA(J1,J2)
2493 sage: (pi_left, pi_right) = J.projections()
2494 sage: J.one().to_vector()
2496 sage: pi_left(J.one()).to_vector()
2498 sage: pi_right(J.one()).to_vector()
2502 (J1
,J2
) = self
.factors()
2505 V_basis
= self
.vector_space().basis()
2506 # Need to specify the dimensions explicitly so that we don't
2507 # wind up with a zero-by-zero matrix when we want e.g. a
2508 # zero-by-two matrix (important for composing things).
2509 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2510 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2511 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2512 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2513 return (pi_left
, pi_right
)
2515 def inclusions(self
):
2517 Return the pair of inclusion maps from our factors into us.
2521 sage: from mjo.eja.eja_algebra import (random_eja,
2522 ....: JordanSpinEJA,
2523 ....: RealSymmetricEJA,
2528 sage: J1 = JordanSpinEJA(3)
2529 sage: J2 = RealSymmetricEJA(2)
2530 sage: J = DirectSumEJA(J1,J2)
2531 sage: (iota_left, iota_right) = J.inclusions()
2532 sage: iota_left(J1.zero()) == J.zero()
2534 sage: iota_right(J2.zero()) == J.zero()
2536 sage: J1.one().to_vector()
2538 sage: iota_left(J1.one()).to_vector()
2540 sage: J2.one().to_vector()
2542 sage: iota_right(J2.one()).to_vector()
2544 sage: J.one().to_vector()
2549 Composing a projection with the corresponding inclusion should
2550 produce the identity map, and mismatching them should produce
2553 sage: set_random_seed()
2554 sage: J1 = random_eja()
2555 sage: J2 = random_eja()
2556 sage: J = DirectSumEJA(J1,J2)
2557 sage: (iota_left, iota_right) = J.inclusions()
2558 sage: (pi_left, pi_right) = J.projections()
2559 sage: pi_left*iota_left == J1.one().operator()
2561 sage: pi_right*iota_right == J2.one().operator()
2563 sage: (pi_left*iota_right).is_zero()
2565 sage: (pi_right*iota_left).is_zero()
2569 (J1
,J2
) = self
.factors()
2572 V_basis
= self
.vector_space().basis()
2573 # Need to specify the dimensions explicitly so that we don't
2574 # wind up with a zero-by-zero matrix when we want e.g. a
2575 # two-by-zero matrix (important for composing things).
2576 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2577 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2578 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2579 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2580 return (iota_left
, iota_right
)
2582 def inner_product(self
, x
, y
):
2584 The standard Cartesian inner-product.
2586 We project ``x`` and ``y`` onto our factors, and add up the
2587 inner-products from the subalgebras.
2592 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2593 ....: QuaternionHermitianEJA,
2598 sage: J1 = HadamardEJA(3,QQ)
2599 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2600 sage: J = DirectSumEJA(J1,J2)
2605 sage: x1.inner_product(x2)
2607 sage: y1.inner_product(y2)
2609 sage: J.one().inner_product(J.one())
2613 (pi_left
, pi_right
) = self
.projections()
2619 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2623 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance