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
.table
import table
28 from sage
.modules
.free_module
import FreeModule
, VectorSpace
29 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
32 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
33 from mjo
.eja
.eja_operator
import FiniteDimensionalEuclideanJordanAlgebraOperator
34 from mjo
.eja
.eja_utils
import _mat2vec
36 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
38 def _coerce_map_from_base_ring(self
):
40 Disable the map from the base ring into the algebra.
42 Performing a nonsense conversion like this automatically
43 is counterpedagogical. The fallback is to try the usual
44 element constructor, which should also fail.
48 sage: from mjo.eja.eja_algebra import random_eja
52 sage: set_random_seed()
53 sage: J = random_eja()
55 Traceback (most recent call last):
57 ValueError: not an element of this algebra
73 sage: from mjo.eja.eja_algebra import (
74 ....: FiniteDimensionalEuclideanJordanAlgebra,
80 By definition, Jordan multiplication commutes::
82 sage: set_random_seed()
83 sage: J = random_eja()
84 sage: x,y = J.random_elements(2)
90 The ``field`` we're given must be real with ``check_field=True``::
92 sage: JordanSpinEJA(2,QQbar)
93 Traceback (most recent call last):
95 ValueError: scalar field is not real
97 The multiplication table must be square with ``check_axioms=True``::
99 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
100 Traceback (most recent call last):
102 ValueError: multiplication table is not square
106 if not field
.is_subring(RR
):
107 # Note: this does return true for the real algebraic
108 # field, the rationals, and any quadratic field where
109 # we've specified a real embedding.
110 raise ValueError("scalar field is not real")
112 # The multiplication table had better be square
115 if not all( len(l
) == n
for l
in mult_table
):
116 raise ValueError("multiplication table is not square")
118 self
._matrix
_basis
= matrix_basis
121 category
= MagmaticAlgebras(field
).FiniteDimensional()
122 category
= category
.WithBasis().Unital()
124 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
129 self
.print_options(bracket
='')
131 # The multiplication table we're given is necessarily in terms
132 # of vectors, because we don't have an algebra yet for
133 # anything to be an element of. However, it's faster in the
134 # long run to have the multiplication table be in terms of
135 # algebra elements. We do this after calling the superclass
136 # constructor so that from_vector() knows what to do.
137 self
._multiplication
_table
= [ [ self
.vector_space().zero()
140 # take advantage of symmetry
143 elt
= self
.from_vector(mult_table
[i
][j
])
144 self
._multiplication
_table
[i
][j
] = elt
145 self
._multiplication
_table
[j
][i
] = elt
148 if not self
._is
_commutative
():
149 raise ValueError("algebra is not commutative")
150 if not self
._is
_jordanian
():
151 raise ValueError("Jordan identity does not hold")
152 if not self
._inner
_product
_is
_associative
():
153 raise ValueError("inner product is not associative")
155 def _element_constructor_(self
, elt
):
157 Construct an element of this algebra from its vector or matrix
160 This gets called only after the parent element _call_ method
161 fails to find a coercion for the argument.
165 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
167 ....: RealSymmetricEJA)
171 The identity in `S^n` is converted to the identity in the EJA::
173 sage: J = RealSymmetricEJA(3)
174 sage: I = matrix.identity(QQ,3)
175 sage: J(I) == J.one()
178 This skew-symmetric matrix can't be represented in the EJA::
180 sage: J = RealSymmetricEJA(3)
181 sage: A = matrix(QQ,3, lambda i,j: i-j)
183 Traceback (most recent call last):
185 ValueError: not an element of this algebra
189 Ensure that we can convert any element of the two non-matrix
190 simple algebras (whose matrix representations are columns)
191 back and forth faithfully::
193 sage: set_random_seed()
194 sage: J = HadamardEJA.random_instance()
195 sage: x = J.random_element()
196 sage: J(x.to_vector().column()) == x
198 sage: J = JordanSpinEJA.random_instance()
199 sage: x = J.random_element()
200 sage: J(x.to_vector().column()) == x
203 msg
= "not an element of this algebra"
205 # The superclass implementation of random_element()
206 # needs to be able to coerce "0" into the algebra.
208 elif elt
in self
.base_ring():
209 # Ensure that no base ring -> algebra coercion is performed
210 # by this method. There's some stupidity in sage that would
211 # otherwise propagate to this method; for example, sage thinks
212 # that the integer 3 belongs to the space of 2-by-2 matrices.
213 raise ValueError(msg
)
215 if elt
not in self
.matrix_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 matrix-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(self
.base_ring(), elt
.nrows()*elt
.ncols())
225 W
= V
.span_of_basis( _mat2vec(s
) for s
in self
.matrix_basis() )
228 coords
= W
.coordinate_vector(_mat2vec(elt
))
229 except ArithmeticError: # vector is not in free module
230 raise ValueError(msg
)
232 return self
.from_vector(coords
)
236 Return a string representation of ``self``.
240 sage: from mjo.eja.eja_algebra import JordanSpinEJA
244 Ensure that it says what we think it says::
246 sage: JordanSpinEJA(2, field=AA)
247 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
248 sage: JordanSpinEJA(3, field=RDF)
249 Euclidean Jordan algebra of dimension 3 over Real Double Field
252 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
253 return fmt
.format(self
.dimension(), self
.base_ring())
255 def product_on_basis(self
, i
, j
):
256 return self
._multiplication
_table
[i
][j
]
258 def _is_commutative(self
):
260 Whether or not this algebra's multiplication table is commutative.
262 This method should of course always return ``True``, unless
263 this algebra was constructed with ``check_axioms=False`` and
264 passed an invalid multiplication table.
266 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
267 for i
in range(self
.dimension())
268 for j
in range(self
.dimension()) )
270 def _is_jordanian(self
):
272 Whether or not this algebra's multiplication table respects the
273 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
275 We only check one arrangement of `x` and `y`, so for a
276 ``True`` result to be truly true, you should also check
277 :meth:`_is_commutative`. This method should of course always
278 return ``True``, unless this algebra was constructed with
279 ``check_axioms=False`` and passed an invalid multiplication table.
281 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
283 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
284 for i
in range(self
.dimension())
285 for j
in range(self
.dimension()) )
287 def _inner_product_is_associative(self
):
289 Return whether or not this algebra's inner product `B` is
290 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
292 This method should of course always return ``True``, unless
293 this algebra was constructed with ``check_axioms=False`` and
294 passed an invalid multiplication table.
297 # Used to check whether or not something is zero in an inexact
298 # ring. This number is sufficient to allow the construction of
299 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
302 for i
in range(self
.dimension()):
303 for j
in range(self
.dimension()):
304 for k
in range(self
.dimension()):
308 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
310 if self
.base_ring().is_exact():
314 if diff
.abs() > epsilon
:
320 def characteristic_polynomial_of(self
):
322 Return the algebra's "characteristic polynomial of" function,
323 which is itself a multivariate polynomial that, when evaluated
324 at the coordinates of some algebra element, returns that
325 element's characteristic polynomial.
327 The resulting polynomial has `n+1` variables, where `n` is the
328 dimension of this algebra. The first `n` variables correspond to
329 the coordinates of an algebra element: when evaluated at the
330 coordinates of an algebra element with respect to a certain
331 basis, the result is a univariate polynomial (in the one
332 remaining variable ``t``), namely the characteristic polynomial
337 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
341 The characteristic polynomial in the spin algebra is given in
342 Alizadeh, Example 11.11::
344 sage: J = JordanSpinEJA(3)
345 sage: p = J.characteristic_polynomial_of(); p
346 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
347 sage: xvec = J.one().to_vector()
351 By definition, the characteristic polynomial is a monic
352 degree-zero polynomial in a rank-zero algebra. Note that
353 Cayley-Hamilton is indeed satisfied since the polynomial
354 ``1`` evaluates to the identity element of the algebra on
357 sage: J = TrivialEJA()
358 sage: J.characteristic_polynomial_of()
365 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
366 a
= self
._charpoly
_coefficients
()
368 # We go to a bit of trouble here to reorder the
369 # indeterminates, so that it's easier to evaluate the
370 # characteristic polynomial at x's coordinates and get back
371 # something in terms of t, which is what we want.
372 S
= PolynomialRing(self
.base_ring(),'t')
376 S
= PolynomialRing(S
, R
.variable_names())
379 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
381 def coordinate_polynomial_ring(self
):
383 The multivariate polynomial ring in which this algebra's
384 :meth:`characteristic_polynomial_of` lives.
388 sage: from mjo.eja.eja_algebra import (HadamardEJA,
389 ....: RealSymmetricEJA)
393 sage: J = HadamardEJA(2)
394 sage: J.coordinate_polynomial_ring()
395 Multivariate Polynomial Ring in X1, X2...
396 sage: J = RealSymmetricEJA(3,QQ)
397 sage: J.coordinate_polynomial_ring()
398 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
401 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
402 return PolynomialRing(self
.base_ring(), var_names
)
404 def inner_product(self
, x
, y
):
406 The inner product associated with this Euclidean Jordan algebra.
408 Defaults to the trace inner product, but can be overridden by
409 subclasses if they are sure that the necessary properties are
414 sage: from mjo.eja.eja_algebra import (random_eja,
416 ....: BilinearFormEJA)
420 Our inner product is "associative," which means the following for
421 a symmetric bilinear form::
423 sage: set_random_seed()
424 sage: J = random_eja()
425 sage: x,y,z = J.random_elements(3)
426 sage: (x*y).inner_product(z) == y.inner_product(x*z)
431 Ensure that this is the usual inner product for the algebras
434 sage: set_random_seed()
435 sage: J = HadamardEJA.random_instance()
436 sage: x,y = J.random_elements(2)
437 sage: actual = x.inner_product(y)
438 sage: expected = x.to_vector().inner_product(y.to_vector())
439 sage: actual == expected
442 Ensure that this is one-half of the trace inner-product in a
443 BilinearFormEJA that isn't just the reals (when ``n`` isn't
444 one). This is in Faraut and Koranyi, and also my "On the
447 sage: set_random_seed()
448 sage: J = BilinearFormEJA.random_instance()
449 sage: n = J.dimension()
450 sage: x = J.random_element()
451 sage: y = J.random_element()
452 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
455 B
= self
._inner
_product
_matrix
456 return (B
*x
.to_vector()).inner_product(y
.to_vector())
459 def is_trivial(self
):
461 Return whether or not this algebra is trivial.
463 A trivial algebra contains only the zero element.
467 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
472 sage: J = ComplexHermitianEJA(3)
478 sage: J = TrivialEJA()
483 return self
.dimension() == 0
486 def multiplication_table(self
):
488 Return a visual representation of this algebra's multiplication
489 table (on basis elements).
493 sage: from mjo.eja.eja_algebra import JordanSpinEJA
497 sage: J = JordanSpinEJA(4)
498 sage: J.multiplication_table()
499 +----++----+----+----+----+
500 | * || e0 | e1 | e2 | e3 |
501 +====++====+====+====+====+
502 | e0 || e0 | e1 | e2 | e3 |
503 +----++----+----+----+----+
504 | e1 || e1 | e0 | 0 | 0 |
505 +----++----+----+----+----+
506 | e2 || e2 | 0 | e0 | 0 |
507 +----++----+----+----+----+
508 | e3 || e3 | 0 | 0 | e0 |
509 +----++----+----+----+----+
512 M
= list(self
._multiplication
_table
) # copy
513 for i
in range(len(M
)):
514 # M had better be "square"
515 M
[i
] = [self
.monomial(i
)] + M
[i
]
516 M
= [["*"] + list(self
.gens())] + M
517 return table(M
, header_row
=True, header_column
=True, frame
=True)
520 def matrix_basis(self
):
522 Return an (often more natural) representation of this algebras
523 basis as an ordered tuple of matrices.
525 Every finite-dimensional Euclidean Jordan Algebra is a, up to
526 Jordan isomorphism, a direct sum of five simple
527 algebras---four of which comprise Hermitian matrices. And the
528 last type of algebra can of course be thought of as `n`-by-`1`
529 column matrices (ambiguusly called column vectors) to avoid
530 special cases. As a result, matrices (and column vectors) are
531 a natural representation format for Euclidean Jordan algebra
534 But, when we construct an algebra from a basis of matrices,
535 those matrix representations are lost in favor of coordinate
536 vectors *with respect to* that basis. We could eventually
537 convert back if we tried hard enough, but having the original
538 representations handy is valuable enough that we simply store
539 them and return them from this method.
541 Why implement this for non-matrix algebras? Avoiding special
542 cases for the :class:`BilinearFormEJA` pays with simplicity in
543 its own right. But mainly, we would like to be able to assume
544 that elements of a :class:`DirectSumEJA` can be displayed
545 nicely, without having to have special classes for direct sums
546 one of whose components was a matrix algebra.
550 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
551 ....: RealSymmetricEJA)
555 sage: J = RealSymmetricEJA(2)
557 Finite family {0: e0, 1: e1, 2: e2}
558 sage: J.matrix_basis()
560 [1 0] [ 0 0.7071067811865475?] [0 0]
561 [0 0], [0.7071067811865475? 0], [0 1]
566 sage: J = JordanSpinEJA(2)
568 Finite family {0: e0, 1: e1}
569 sage: J.matrix_basis()
575 if self
._matrix
_basis
is None:
576 M
= self
.matrix_space()
577 return tuple( M(b
.to_vector()) for b
in self
.basis() )
579 return self
._matrix
_basis
582 def matrix_space(self
):
584 Return the matrix space in which this algebra's elements live, if
585 we think of them as matrices (including column vectors of the
588 Generally this will be an `n`-by-`1` column-vector space,
589 except when the algebra is trivial. There it's `n`-by-`n`
590 (where `n` is zero), to ensure that two elements of the matrix
591 space (empty matrices) can be multiplied.
593 Matrix algebras override this with something more useful.
595 if self
.is_trivial():
596 return MatrixSpace(self
.base_ring(), 0)
597 elif self
._matrix
_basis
is None or len(self
._matrix
_basis
) == 0:
598 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
600 return self
._matrix
_basis
[0].matrix_space()
606 Return the unit element of this algebra.
610 sage: from mjo.eja.eja_algebra import (HadamardEJA,
615 sage: J = HadamardEJA(5)
617 e0 + e1 + e2 + e3 + e4
621 The identity element acts like the identity::
623 sage: set_random_seed()
624 sage: J = random_eja()
625 sage: x = J.random_element()
626 sage: J.one()*x == x and x*J.one() == x
629 The matrix of the unit element's operator is the identity::
631 sage: set_random_seed()
632 sage: J = random_eja()
633 sage: actual = J.one().operator().matrix()
634 sage: expected = matrix.identity(J.base_ring(), J.dimension())
635 sage: actual == expected
638 Ensure that the cached unit element (often precomputed by
639 hand) agrees with the computed one::
641 sage: set_random_seed()
642 sage: J = random_eja()
643 sage: cached = J.one()
644 sage: J.one.clear_cache()
645 sage: J.one() == cached
649 # We can brute-force compute the matrices of the operators
650 # that correspond to the basis elements of this algebra.
651 # If some linear combination of those basis elements is the
652 # algebra identity, then the same linear combination of
653 # their matrices has to be the identity matrix.
655 # Of course, matrices aren't vectors in sage, so we have to
656 # appeal to the "long vectors" isometry.
657 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
659 # Now we use basic linear algebra to find the coefficients,
660 # of the matrices-as-vectors-linear-combination, which should
661 # work for the original algebra basis too.
662 A
= matrix(self
.base_ring(), oper_vecs
)
664 # We used the isometry on the left-hand side already, but we
665 # still need to do it for the right-hand side. Recall that we
666 # wanted something that summed to the identity matrix.
667 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
669 # Now if there's an identity element in the algebra, this
670 # should work. We solve on the left to avoid having to
671 # transpose the matrix "A".
672 return self
.from_vector(A
.solve_left(b
))
675 def peirce_decomposition(self
, c
):
677 The Peirce decomposition of this algebra relative to the
680 In the future, this can be extended to a complete system of
681 orthogonal idempotents.
685 - ``c`` -- an idempotent of this algebra.
689 A triple (J0, J5, J1) containing two subalgebras and one subspace
692 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
693 corresponding to the eigenvalue zero.
695 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
696 corresponding to the eigenvalue one-half.
698 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
699 corresponding to the eigenvalue one.
701 These are the only possible eigenspaces for that operator, and this
702 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
703 orthogonal, and are subalgebras of this algebra with the appropriate
708 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
712 The canonical example comes from the symmetric matrices, which
713 decompose into diagonal and off-diagonal parts::
715 sage: J = RealSymmetricEJA(3)
716 sage: C = matrix(QQ, [ [1,0,0],
720 sage: J0,J5,J1 = J.peirce_decomposition(c)
722 Euclidean Jordan algebra of dimension 1...
724 Vector space of degree 6 and dimension 2...
726 Euclidean Jordan algebra of dimension 3...
727 sage: J0.one().to_matrix()
731 sage: orig_df = AA.options.display_format
732 sage: AA.options.display_format = 'radical'
733 sage: J.from_vector(J5.basis()[0]).to_matrix()
737 sage: J.from_vector(J5.basis()[1]).to_matrix()
741 sage: AA.options.display_format = orig_df
742 sage: J1.one().to_matrix()
749 Every algebra decomposes trivially with respect to its identity
752 sage: set_random_seed()
753 sage: J = random_eja()
754 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
755 sage: J0.dimension() == 0 and J5.dimension() == 0
757 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
760 The decomposition is into eigenspaces, and its components are
761 therefore necessarily orthogonal. Moreover, the identity
762 elements in the two subalgebras are the projections onto their
763 respective subspaces of the superalgebra's identity element::
765 sage: set_random_seed()
766 sage: J = random_eja()
767 sage: x = J.random_element()
768 sage: if not J.is_trivial():
769 ....: while x.is_nilpotent():
770 ....: x = J.random_element()
771 sage: c = x.subalgebra_idempotent()
772 sage: J0,J5,J1 = J.peirce_decomposition(c)
774 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
775 ....: w = w.superalgebra_element()
776 ....: y = J.from_vector(y)
777 ....: z = z.superalgebra_element()
778 ....: ipsum += w.inner_product(y).abs()
779 ....: ipsum += w.inner_product(z).abs()
780 ....: ipsum += y.inner_product(z).abs()
783 sage: J1(c) == J1.one()
785 sage: J0(J.one() - c) == J0.one()
789 if not c
.is_idempotent():
790 raise ValueError("element is not idempotent: %s" % c
)
792 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
794 # Default these to what they should be if they turn out to be
795 # trivial, because eigenspaces_left() won't return eigenvalues
796 # corresponding to trivial spaces (e.g. it returns only the
797 # eigenspace corresponding to lambda=1 if you take the
798 # decomposition relative to the identity element).
799 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
800 J0
= trivial
# eigenvalue zero
801 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
802 J1
= trivial
# eigenvalue one
804 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
805 if eigval
== ~
(self
.base_ring()(2)):
808 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
809 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
817 raise ValueError("unexpected eigenvalue: %s" % eigval
)
822 def random_element(self
, thorough
=False):
824 Return a random element of this algebra.
826 Our algebra superclass method only returns a linear
827 combination of at most two basis elements. We instead
828 want the vector space "random element" method that
829 returns a more diverse selection.
833 - ``thorough`` -- (boolean; default False) whether or not we
834 should generate irrational coefficients for the random
835 element when our base ring is irrational; this slows the
836 algebra operations to a crawl, but any truly random method
840 # For a general base ring... maybe we can trust this to do the
841 # right thing? Unlikely, but.
842 V
= self
.vector_space()
843 v
= V
.random_element()
845 if self
.base_ring() is AA
:
846 # The "random element" method of the algebraic reals is
847 # stupid at the moment, and only returns integers between
848 # -2 and 2, inclusive:
850 # https://trac.sagemath.org/ticket/30875
852 # Instead, we implement our own "random vector" method,
853 # and then coerce that into the algebra. We use the vector
854 # space degree here instead of the dimension because a
855 # subalgebra could (for example) be spanned by only two
856 # vectors, each with five coordinates. We need to
857 # generate all five coordinates.
859 v
*= QQbar
.random_element().real()
861 v
*= QQ
.random_element()
863 return self
.from_vector(V
.coordinate_vector(v
))
865 def random_elements(self
, count
, thorough
=False):
867 Return ``count`` random elements as a tuple.
871 - ``thorough`` -- (boolean; default False) whether or not we
872 should generate irrational coefficients for the random
873 elements when our base ring is irrational; this slows the
874 algebra operations to a crawl, but any truly random method
879 sage: from mjo.eja.eja_algebra import JordanSpinEJA
883 sage: J = JordanSpinEJA(3)
884 sage: x,y,z = J.random_elements(3)
885 sage: all( [ x in J, y in J, z in J ])
887 sage: len( J.random_elements(10) ) == 10
891 return tuple( self
.random_element(thorough
)
892 for idx
in range(count
) )
896 def _charpoly_coefficients(self
):
898 The `r` polynomial coefficients of the "characteristic polynomial
902 R
= self
.coordinate_polynomial_ring()
904 F
= R
.fraction_field()
907 # From a result in my book, these are the entries of the
908 # basis representation of L_x.
909 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
912 L_x
= matrix(F
, n
, n
, L_x_i_j
)
915 if self
.rank
.is_in_cache():
917 # There's no need to pad the system with redundant
918 # columns if we *know* they'll be redundant.
921 # Compute an extra power in case the rank is equal to
922 # the dimension (otherwise, we would stop at x^(r-1)).
923 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
924 for k
in range(n
+1) ]
925 A
= matrix
.column(F
, x_powers
[:n
])
926 AE
= A
.extended_echelon_form()
933 # The theory says that only the first "r" coefficients are
934 # nonzero, and they actually live in the original polynomial
935 # ring and not the fraction field. We negate them because
936 # in the actual characteristic polynomial, they get moved
937 # to the other side where x^r lives.
938 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
943 Return the rank of this EJA.
945 This is a cached method because we know the rank a priori for
946 all of the algebras we can construct. Thus we can avoid the
947 expensive ``_charpoly_coefficients()`` call unless we truly
948 need to compute the whole characteristic polynomial.
952 sage: from mjo.eja.eja_algebra import (HadamardEJA,
954 ....: RealSymmetricEJA,
955 ....: ComplexHermitianEJA,
956 ....: QuaternionHermitianEJA,
961 The rank of the Jordan spin algebra is always two::
963 sage: JordanSpinEJA(2).rank()
965 sage: JordanSpinEJA(3).rank()
967 sage: JordanSpinEJA(4).rank()
970 The rank of the `n`-by-`n` Hermitian real, complex, or
971 quaternion matrices is `n`::
973 sage: RealSymmetricEJA(4).rank()
975 sage: ComplexHermitianEJA(3).rank()
977 sage: QuaternionHermitianEJA(2).rank()
982 Ensure that every EJA that we know how to construct has a
983 positive integer rank, unless the algebra is trivial in
984 which case its rank will be zero::
986 sage: set_random_seed()
987 sage: J = random_eja()
991 sage: r > 0 or (r == 0 and J.is_trivial())
994 Ensure that computing the rank actually works, since the ranks
995 of all simple algebras are known and will be cached by default::
997 sage: set_random_seed() # long time
998 sage: J = random_eja() # long time
999 sage: caches = J.rank() # long time
1000 sage: J.rank.clear_cache() # long time
1001 sage: J.rank() == cached # long time
1005 return len(self
._charpoly
_coefficients
())
1008 def vector_space(self
):
1010 Return the vector space that underlies this algebra.
1014 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1018 sage: J = RealSymmetricEJA(2)
1019 sage: J.vector_space()
1020 Vector space of dimension 3 over...
1023 return self
.zero().to_vector().parent().ambient_vector_space()
1026 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1029 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1031 Algebras whose basis consists of vectors with rational
1032 entries. Equivalently, algebras whose multiplication tables
1033 contain only rational coefficients.
1035 When an EJA has a basis that can be made rational, we can speed up
1036 the computation of its characteristic polynomial by doing it over
1037 ``QQ``. All of the named EJA constructors that we provide fall
1041 def _charpoly_coefficients(self
):
1043 Override the parent method with something that tries to compute
1044 over a faster (non-extension) field.
1048 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1052 The base ring of the resulting polynomial coefficients is what
1053 it should be, and not the rationals (unless the algebra was
1054 already over the rationals)::
1056 sage: J = JordanSpinEJA(3)
1057 sage: J._charpoly_coefficients()
1058 (X1^2 - X2^2 - X3^2, -2*X1)
1059 sage: a0 = J._charpoly_coefficients()[0]
1061 Algebraic Real Field
1062 sage: a0.base_ring()
1063 Algebraic Real Field
1066 if self
.base_ring() is QQ
:
1067 # There's no need to construct *another* algebra over the
1068 # rationals if this one is already over the rationals.
1069 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1070 return superclass
._charpoly
_coefficients
()
1073 tuple(map(lambda x
: x
.to_vector(), ls
))
1074 for ls
in self
._multiplication
_table
1077 # Do the computation over the rationals. The answer will be
1078 # the same, because our basis coordinates are (essentially)
1080 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1084 a
= J
._charpoly
_coefficients
()
1085 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1088 class ConcreteEuclideanJordanAlgebra
:
1090 A class for the Euclidean Jordan algebras that we know by name.
1092 These are the Jordan algebras whose basis, multiplication table,
1093 rank, and so on are known a priori. More to the point, they are
1094 the Euclidean Jordan algebras for which we are able to conjure up
1095 a "random instance."
1099 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1103 Our basis is normalized with respect to the algebra's inner
1104 product, unless we specify otherwise::
1106 sage: set_random_seed()
1107 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1108 sage: all( b.norm() == 1 for b in J.gens() )
1111 Since our basis is orthonormal with respect to the algebra's inner
1112 product, and since we know that this algebra is an EJA, any
1113 left-multiplication operator's matrix will be symmetric because
1114 natural->EJA basis representation is an isometry and within the
1115 EJA the operator is self-adjoint by the Jordan axiom::
1117 sage: set_random_seed()
1118 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1119 sage: x = J.random_element()
1120 sage: x.operator().is_self_adjoint()
1125 def _max_random_instance_size():
1127 Return an integer "size" that is an upper bound on the size of
1128 this algebra when it is used in a random test
1129 case. Unfortunately, the term "size" is ambiguous -- when
1130 dealing with `R^n` under either the Hadamard or Jordan spin
1131 product, the "size" refers to the dimension `n`. When dealing
1132 with a matrix algebra (real symmetric or complex/quaternion
1133 Hermitian), it refers to the size of the matrix, which is far
1134 less than the dimension of the underlying vector space.
1136 This method must be implemented in each subclass.
1138 raise NotImplementedError
1141 def random_instance(cls
, field
=AA
, **kwargs
):
1143 Return a random instance of this type of algebra.
1145 This method should be implemented in each subclass.
1147 from sage
.misc
.prandom
import choice
1148 eja_class
= choice(cls
.__subclasses
__())
1149 return eja_class
.random_instance(field
)
1152 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1154 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1156 Compared to the superclass constructor, we take a basis instead of
1157 a multiplication table because the latter can be computed in terms
1158 of the former when the product is known (like it is here).
1160 # Used in this class's fast _charpoly_coefficients() override.
1161 self
._basis
_normalizers
= None
1163 # We're going to loop through this a few times, so now's a good
1164 # time to ensure that it isn't a generator expression.
1165 basis
= tuple(basis
)
1167 algebra_dim
= len(basis
)
1168 degree
= 0 # size of the matrices
1170 degree
= basis
[0].nrows()
1172 if algebra_dim
> 1 and normalize_basis
:
1173 # We'll need sqrt(2) to normalize the basis, and this
1174 # winds up in the multiplication table, so the whole
1175 # algebra needs to be over the field extension.
1176 R
= PolynomialRing(field
, 'z')
1179 if p
.is_irreducible():
1180 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1181 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1182 self
._basis
_normalizers
= tuple(
1183 ~
(self
.matrix_inner_product(s
,s
).sqrt()) for s
in basis
)
1184 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1186 # Now compute the multiplication and inner product tables.
1187 # We have to do this *after* normalizing the basis, because
1188 # scaling affects the answers.
1189 V
= VectorSpace(field
, degree
**2)
1190 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1191 mult_table
= [[W
.zero() for j
in range(algebra_dim
)]
1192 for i
in range(algebra_dim
)]
1193 ip_table
= [[W
.zero() for j
in range(algebra_dim
)]
1194 for i
in range(algebra_dim
)]
1195 for i
in range(algebra_dim
):
1196 for j
in range(algebra_dim
):
1197 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1198 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1201 # HACK: ignore the error here if we don't need the
1202 # inner product (as is the case when we construct
1203 # a dummy QQ-algebra for fast charpoly coefficients.
1204 ip_table
[i
][j
] = self
.matrix_inner_product(basis
[i
],
1211 self
._inner
_product
_matrix
= matrix(field
,ip_table
)
1215 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1220 if algebra_dim
== 0:
1221 self
.one
.set_cache(self
.zero())
1223 n
= basis
[0].nrows()
1224 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1225 # details of this algebra.
1226 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1230 def _charpoly_coefficients(self
):
1232 Override the parent method with something that tries to compute
1233 over a faster (non-extension) field.
1235 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1236 # We didn't normalize, or the basis we started with had
1237 # entries in a nice field already. Just compute the thing.
1238 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1240 basis
= ( (b
/n
) for (b
,n
) in zip(self
.matrix_basis(),
1241 self
._basis
_normalizers
) )
1243 # Do this over the rationals and convert back at the end.
1244 # Only works because we know the entries of the basis are
1245 # integers. The argument ``check_axioms=False`` is required
1246 # because the trace inner-product method for this
1247 # class is a stub and can't actually be checked.
1248 J
= MatrixEuclideanJordanAlgebra(QQ
,
1250 normalize_basis
=False,
1253 a
= J
._charpoly
_coefficients
()
1255 # Unfortunately, changing the basis does change the
1256 # coefficients of the characteristic polynomial, but since
1257 # these are really the coefficients of the "characteristic
1258 # polynomial of" function, everything is still nice and
1259 # unevaluated. It's therefore "obvious" how scaling the
1260 # basis affects the coordinate variables X1, X2, et
1261 # cetera. Scaling the first basis vector up by "n" adds a
1262 # factor of 1/n into every "X1" term, for example. So here
1263 # we simply undo the basis_normalizer scaling that we
1264 # performed earlier.
1266 # The a[0] access here is safe because trivial algebras
1267 # won't have any basis normalizers and therefore won't
1268 # make it to this "else" branch.
1269 XS
= a
[0].parent().gens()
1270 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1271 for i
in range(len(XS
)) }
1272 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1278 Embed the matrix ``M`` into a space of real matrices.
1280 The matrix ``M`` can have entries in any field at the moment:
1281 the real numbers, complex numbers, or quaternions. And although
1282 they are not a field, we can probably support octonions at some
1283 point, too. This function returns a real matrix that "acts like"
1284 the original with respect to matrix multiplication; i.e.
1286 real_embed(M*N) = real_embed(M)*real_embed(N)
1289 raise NotImplementedError
1293 def real_unembed(M
):
1295 The inverse of :meth:`real_embed`.
1297 raise NotImplementedError
1300 def matrix_inner_product(cls
,X
,Y
):
1301 Xu
= cls
.real_unembed(X
)
1302 Yu
= cls
.real_unembed(Y
)
1303 tr
= (Xu
*Yu
).trace()
1306 # Works in QQ, AA, RDF, et cetera.
1308 except AttributeError:
1309 # A quaternion doesn't have a real() method, but does
1310 # have coefficient_tuple() method that returns the
1311 # coefficients of 1, i, j, and k -- in that order.
1312 return tr
.coefficient_tuple()[0]
1315 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1319 The identity function, for embedding real matrices into real
1325 def real_unembed(M
):
1327 The identity function, for unembedding real matrices from real
1333 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1334 ConcreteEuclideanJordanAlgebra
):
1336 The rank-n simple EJA consisting of real symmetric n-by-n
1337 matrices, the usual symmetric Jordan product, and the trace inner
1338 product. It has dimension `(n^2 + n)/2` over the reals.
1342 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1346 sage: J = RealSymmetricEJA(2)
1347 sage: e0, e1, e2 = J.gens()
1355 In theory, our "field" can be any subfield of the reals::
1357 sage: RealSymmetricEJA(2, RDF)
1358 Euclidean Jordan algebra of dimension 3 over Real Double Field
1359 sage: RealSymmetricEJA(2, RR)
1360 Euclidean Jordan algebra of dimension 3 over Real Field with
1361 53 bits of precision
1365 The dimension of this algebra is `(n^2 + n) / 2`::
1367 sage: set_random_seed()
1368 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1369 sage: n = ZZ.random_element(1, n_max)
1370 sage: J = RealSymmetricEJA(n)
1371 sage: J.dimension() == (n^2 + n)/2
1374 The Jordan multiplication is what we think it is::
1376 sage: set_random_seed()
1377 sage: J = RealSymmetricEJA.random_instance()
1378 sage: x,y = J.random_elements(2)
1379 sage: actual = (x*y).to_matrix()
1380 sage: X = x.to_matrix()
1381 sage: Y = y.to_matrix()
1382 sage: expected = (X*Y + Y*X)/2
1383 sage: actual == expected
1385 sage: J(expected) == x*y
1388 We can change the generator prefix::
1390 sage: RealSymmetricEJA(3, prefix='q').gens()
1391 (q0, q1, q2, q3, q4, q5)
1393 We can construct the (trivial) algebra of rank zero::
1395 sage: RealSymmetricEJA(0)
1396 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1400 def _denormalized_basis(cls
, n
, field
):
1402 Return a basis for the space of real symmetric n-by-n matrices.
1406 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1410 sage: set_random_seed()
1411 sage: n = ZZ.random_element(1,5)
1412 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1413 sage: all( M.is_symmetric() for M in B)
1417 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1421 for j
in range(i
+1):
1422 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1426 Sij
= Eij
+ Eij
.transpose()
1432 def _max_random_instance_size():
1433 return 4 # Dimension 10
1436 def random_instance(cls
, field
=AA
, **kwargs
):
1438 Return a random instance of this type of algebra.
1440 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1441 return cls(n
, field
, **kwargs
)
1443 def __init__(self
, n
, field
=AA
, **kwargs
):
1444 basis
= self
._denormalized
_basis
(n
, field
)
1445 super(RealSymmetricEJA
, self
).__init
__(field
,
1449 self
.rank
.set_cache(n
)
1452 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1456 Embed the n-by-n complex matrix ``M`` into the space of real
1457 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1458 bi` to the block matrix ``[[a,b],[-b,a]]``.
1462 sage: from mjo.eja.eja_algebra import \
1463 ....: ComplexMatrixEuclideanJordanAlgebra
1467 sage: F = QuadraticField(-1, 'I')
1468 sage: x1 = F(4 - 2*i)
1469 sage: x2 = F(1 + 2*i)
1472 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1473 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1482 Embedding is a homomorphism (isomorphism, in fact)::
1484 sage: set_random_seed()
1485 sage: n = ZZ.random_element(3)
1486 sage: F = QuadraticField(-1, 'I')
1487 sage: X = random_matrix(F, n)
1488 sage: Y = random_matrix(F, n)
1489 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1490 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1491 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1498 raise ValueError("the matrix 'M' must be square")
1500 # We don't need any adjoined elements...
1501 field
= M
.base_ring().base_ring()
1505 a
= z
.list()[0] # real part, I guess
1506 b
= z
.list()[1] # imag part, I guess
1507 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1509 return matrix
.block(field
, n
, blocks
)
1513 def real_unembed(M
):
1515 The inverse of _embed_complex_matrix().
1519 sage: from mjo.eja.eja_algebra import \
1520 ....: ComplexMatrixEuclideanJordanAlgebra
1524 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1525 ....: [-2, 1, -4, 3],
1526 ....: [ 9, 10, 11, 12],
1527 ....: [-10, 9, -12, 11] ])
1528 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1530 [ 10*I + 9 12*I + 11]
1534 Unembedding is the inverse of embedding::
1536 sage: set_random_seed()
1537 sage: F = QuadraticField(-1, 'I')
1538 sage: M = random_matrix(F, 3)
1539 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1540 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1546 raise ValueError("the matrix 'M' must be square")
1547 if not n
.mod(2).is_zero():
1548 raise ValueError("the matrix 'M' must be a complex embedding")
1550 # If "M" was normalized, its base ring might have roots
1551 # adjoined and they can stick around after unembedding.
1552 field
= M
.base_ring()
1553 R
= PolynomialRing(field
, 'z')
1556 # Sage doesn't know how to embed AA into QQbar, i.e. how
1557 # to adjoin sqrt(-1) to AA.
1560 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1563 # Go top-left to bottom-right (reading order), converting every
1564 # 2-by-2 block we see to a single complex element.
1566 for k
in range(n
/2):
1567 for j
in range(n
/2):
1568 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1569 if submat
[0,0] != submat
[1,1]:
1570 raise ValueError('bad on-diagonal submatrix')
1571 if submat
[0,1] != -submat
[1,0]:
1572 raise ValueError('bad off-diagonal submatrix')
1573 z
= submat
[0,0] + submat
[0,1]*i
1576 return matrix(F
, n
/2, elements
)
1580 def matrix_inner_product(cls
,X
,Y
):
1582 Compute a matrix inner product in this algebra directly from
1587 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1591 This gives the same answer as the slow, default method implemented
1592 in :class:`MatrixEuclideanJordanAlgebra`::
1594 sage: set_random_seed()
1595 sage: J = ComplexHermitianEJA.random_instance()
1596 sage: x,y = J.random_elements(2)
1597 sage: Xe = x.to_matrix()
1598 sage: Ye = y.to_matrix()
1599 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1600 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1601 sage: expected = (X*Y).trace().real()
1602 sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
1603 sage: actual == expected
1607 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/2
1610 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1611 ConcreteEuclideanJordanAlgebra
):
1613 The rank-n simple EJA consisting of complex Hermitian n-by-n
1614 matrices over the real numbers, the usual symmetric Jordan product,
1615 and the real-part-of-trace inner product. It has dimension `n^2` over
1620 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1624 In theory, our "field" can be any subfield of the reals::
1626 sage: ComplexHermitianEJA(2, RDF)
1627 Euclidean Jordan algebra of dimension 4 over Real Double Field
1628 sage: ComplexHermitianEJA(2, RR)
1629 Euclidean Jordan algebra of dimension 4 over Real Field with
1630 53 bits of precision
1634 The dimension of this algebra is `n^2`::
1636 sage: set_random_seed()
1637 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1638 sage: n = ZZ.random_element(1, n_max)
1639 sage: J = ComplexHermitianEJA(n)
1640 sage: J.dimension() == n^2
1643 The Jordan multiplication is what we think it is::
1645 sage: set_random_seed()
1646 sage: J = ComplexHermitianEJA.random_instance()
1647 sage: x,y = J.random_elements(2)
1648 sage: actual = (x*y).to_matrix()
1649 sage: X = x.to_matrix()
1650 sage: Y = y.to_matrix()
1651 sage: expected = (X*Y + Y*X)/2
1652 sage: actual == expected
1654 sage: J(expected) == x*y
1657 We can change the generator prefix::
1659 sage: ComplexHermitianEJA(2, prefix='z').gens()
1662 We can construct the (trivial) algebra of rank zero::
1664 sage: ComplexHermitianEJA(0)
1665 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1670 def _denormalized_basis(cls
, n
, field
):
1672 Returns a basis for the space of complex Hermitian n-by-n matrices.
1674 Why do we embed these? Basically, because all of numerical linear
1675 algebra assumes that you're working with vectors consisting of `n`
1676 entries from a field and scalars from the same field. There's no way
1677 to tell SageMath that (for example) the vectors contain complex
1678 numbers, while the scalar field is real.
1682 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1686 sage: set_random_seed()
1687 sage: n = ZZ.random_element(1,5)
1688 sage: field = QuadraticField(2, 'sqrt2')
1689 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1690 sage: all( M.is_symmetric() for M in B)
1694 R
= PolynomialRing(field
, 'z')
1696 F
= field
.extension(z
**2 + 1, 'I')
1699 # This is like the symmetric case, but we need to be careful:
1701 # * We want conjugate-symmetry, not just symmetry.
1702 # * The diagonal will (as a result) be real.
1706 for j
in range(i
+1):
1707 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1709 Sij
= cls
.real_embed(Eij
)
1712 # The second one has a minus because it's conjugated.
1713 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1715 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1718 # Since we embedded these, we can drop back to the "field" that we
1719 # started with instead of the complex extension "F".
1720 return ( s
.change_ring(field
) for s
in S
)
1723 def __init__(self
, n
, field
=AA
, **kwargs
):
1724 basis
= self
._denormalized
_basis
(n
,field
)
1725 super(ComplexHermitianEJA
,self
).__init
__(field
,
1729 self
.rank
.set_cache(n
)
1732 def _max_random_instance_size():
1733 return 3 # Dimension 9
1736 def random_instance(cls
, field
=AA
, **kwargs
):
1738 Return a random instance of this type of algebra.
1740 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1741 return cls(n
, field
, **kwargs
)
1743 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1747 Embed the n-by-n quaternion matrix ``M`` into the space of real
1748 matrices of size 4n-by-4n by first sending each quaternion entry `z
1749 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1750 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1755 sage: from mjo.eja.eja_algebra import \
1756 ....: QuaternionMatrixEuclideanJordanAlgebra
1760 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1761 sage: i,j,k = Q.gens()
1762 sage: x = 1 + 2*i + 3*j + 4*k
1763 sage: M = matrix(Q, 1, [[x]])
1764 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1770 Embedding is a homomorphism (isomorphism, in fact)::
1772 sage: set_random_seed()
1773 sage: n = ZZ.random_element(2)
1774 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1775 sage: X = random_matrix(Q, n)
1776 sage: Y = random_matrix(Q, n)
1777 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1778 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1779 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1784 quaternions
= M
.base_ring()
1787 raise ValueError("the matrix 'M' must be square")
1789 F
= QuadraticField(-1, 'I')
1794 t
= z
.coefficient_tuple()
1799 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1800 [-c
+ d
*i
, a
- b
*i
]])
1801 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1802 blocks
.append(realM
)
1804 # We should have real entries by now, so use the realest field
1805 # we've got for the return value.
1806 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1811 def real_unembed(M
):
1813 The inverse of _embed_quaternion_matrix().
1817 sage: from mjo.eja.eja_algebra import \
1818 ....: QuaternionMatrixEuclideanJordanAlgebra
1822 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1823 ....: [-2, 1, -4, 3],
1824 ....: [-3, 4, 1, -2],
1825 ....: [-4, -3, 2, 1]])
1826 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1827 [1 + 2*i + 3*j + 4*k]
1831 Unembedding is the inverse of embedding::
1833 sage: set_random_seed()
1834 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1835 sage: M = random_matrix(Q, 3)
1836 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1837 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1843 raise ValueError("the matrix 'M' must be square")
1844 if not n
.mod(4).is_zero():
1845 raise ValueError("the matrix 'M' must be a quaternion embedding")
1847 # Use the base ring of the matrix to ensure that its entries can be
1848 # multiplied by elements of the quaternion algebra.
1849 field
= M
.base_ring()
1850 Q
= QuaternionAlgebra(field
,-1,-1)
1853 # Go top-left to bottom-right (reading order), converting every
1854 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1857 for l
in range(n
/4):
1858 for m
in range(n
/4):
1859 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1860 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1861 if submat
[0,0] != submat
[1,1].conjugate():
1862 raise ValueError('bad on-diagonal submatrix')
1863 if submat
[0,1] != -submat
[1,0].conjugate():
1864 raise ValueError('bad off-diagonal submatrix')
1865 z
= submat
[0,0].real()
1866 z
+= submat
[0,0].imag()*i
1867 z
+= submat
[0,1].real()*j
1868 z
+= submat
[0,1].imag()*k
1871 return matrix(Q
, n
/4, elements
)
1875 def matrix_inner_product(cls
,X
,Y
):
1877 Compute a matrix inner product in this algebra directly from
1882 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1886 This gives the same answer as the slow, default method implemented
1887 in :class:`MatrixEuclideanJordanAlgebra`::
1889 sage: set_random_seed()
1890 sage: J = QuaternionHermitianEJA.random_instance()
1891 sage: x,y = J.random_elements(2)
1892 sage: Xe = x.to_matrix()
1893 sage: Ye = y.to_matrix()
1894 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1895 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1896 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1897 sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
1898 sage: actual == expected
1902 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/4
1905 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1906 ConcreteEuclideanJordanAlgebra
):
1908 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1909 matrices, the usual symmetric Jordan product, and the
1910 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1915 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1919 In theory, our "field" can be any subfield of the reals::
1921 sage: QuaternionHermitianEJA(2, RDF)
1922 Euclidean Jordan algebra of dimension 6 over Real Double Field
1923 sage: QuaternionHermitianEJA(2, RR)
1924 Euclidean Jordan algebra of dimension 6 over Real Field with
1925 53 bits of precision
1929 The dimension of this algebra is `2*n^2 - n`::
1931 sage: set_random_seed()
1932 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
1933 sage: n = ZZ.random_element(1, n_max)
1934 sage: J = QuaternionHermitianEJA(n)
1935 sage: J.dimension() == 2*(n^2) - n
1938 The Jordan multiplication is what we think it is::
1940 sage: set_random_seed()
1941 sage: J = QuaternionHermitianEJA.random_instance()
1942 sage: x,y = J.random_elements(2)
1943 sage: actual = (x*y).to_matrix()
1944 sage: X = x.to_matrix()
1945 sage: Y = y.to_matrix()
1946 sage: expected = (X*Y + Y*X)/2
1947 sage: actual == expected
1949 sage: J(expected) == x*y
1952 We can change the generator prefix::
1954 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1955 (a0, a1, a2, a3, a4, a5)
1957 We can construct the (trivial) algebra of rank zero::
1959 sage: QuaternionHermitianEJA(0)
1960 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1964 def _denormalized_basis(cls
, n
, field
):
1966 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1968 Why do we embed these? Basically, because all of numerical
1969 linear algebra assumes that you're working with vectors consisting
1970 of `n` entries from a field and scalars from the same field. There's
1971 no way to tell SageMath that (for example) the vectors contain
1972 complex numbers, while the scalar field is real.
1976 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1980 sage: set_random_seed()
1981 sage: n = ZZ.random_element(1,5)
1982 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1983 sage: all( M.is_symmetric() for M in B )
1987 Q
= QuaternionAlgebra(QQ
,-1,-1)
1990 # This is like the symmetric case, but we need to be careful:
1992 # * We want conjugate-symmetry, not just symmetry.
1993 # * The diagonal will (as a result) be real.
1997 for j
in range(i
+1):
1998 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2000 Sij
= cls
.real_embed(Eij
)
2003 # The second, third, and fourth ones have a minus
2004 # because they're conjugated.
2005 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2007 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2009 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2011 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2014 # Since we embedded these, we can drop back to the "field" that we
2015 # started with instead of the quaternion algebra "Q".
2016 return ( s
.change_ring(field
) for s
in S
)
2019 def __init__(self
, n
, field
=AA
, **kwargs
):
2020 basis
= self
._denormalized
_basis
(n
,field
)
2021 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2025 self
.rank
.set_cache(n
)
2028 def _max_random_instance_size():
2030 The maximum rank of a random QuaternionHermitianEJA.
2032 return 2 # Dimension 6
2035 def random_instance(cls
, field
=AA
, **kwargs
):
2037 Return a random instance of this type of algebra.
2039 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2040 return cls(n
, field
, **kwargs
)
2043 class HadamardEJA(RationalBasisEuclideanJordanAlgebra
,
2044 ConcreteEuclideanJordanAlgebra
):
2046 Return the Euclidean Jordan Algebra corresponding to the set
2047 `R^n` under the Hadamard product.
2049 Note: this is nothing more than the Cartesian product of ``n``
2050 copies of the spin algebra. Once Cartesian product algebras
2051 are implemented, this can go.
2055 sage: from mjo.eja.eja_algebra import HadamardEJA
2059 This multiplication table can be verified by hand::
2061 sage: J = HadamardEJA(3)
2062 sage: e0,e1,e2 = J.gens()
2078 We can change the generator prefix::
2080 sage: HadamardEJA(3, prefix='r').gens()
2084 def __init__(self
, n
, field
=AA
, **kwargs
):
2085 V
= VectorSpace(field
, n
)
2086 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
2089 # Inner products are real numbers and not algebra
2090 # elements, so once we turn the algebra element
2091 # into a vector in inner_product(), we never go
2092 # back. As a result -- contrary to what we do with
2093 # self._multiplication_table -- we store the inner
2094 # product table as a plain old matrix and not as
2095 # an algebra operator.
2096 ip_table
= matrix
.identity(field
,n
)
2097 self
._inner
_product
_matrix
= ip_table
2099 super(HadamardEJA
, self
).__init
__(field
,
2103 self
.rank
.set_cache(n
)
2106 self
.one
.set_cache( self
.zero() )
2108 self
.one
.set_cache( sum(self
.gens()) )
2111 def _max_random_instance_size():
2113 The maximum dimension of a random HadamardEJA.
2118 def random_instance(cls
, field
=AA
, **kwargs
):
2120 Return a random instance of this type of algebra.
2122 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2123 return cls(n
, field
, **kwargs
)
2126 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
,
2127 ConcreteEuclideanJordanAlgebra
):
2129 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2130 with the half-trace inner product and jordan product ``x*y =
2131 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2132 a symmetric positive-definite "bilinear form" matrix. Its
2133 dimension is the size of `B`, and it has rank two in dimensions
2134 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2135 the identity matrix of order ``n``.
2137 We insist that the one-by-one upper-left identity block of `B` be
2138 passed in as well so that we can be passed a matrix of size zero
2139 to construct a trivial algebra.
2143 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2144 ....: JordanSpinEJA)
2148 When no bilinear form is specified, the identity matrix is used,
2149 and the resulting algebra is the Jordan spin algebra::
2151 sage: B = matrix.identity(AA,3)
2152 sage: J0 = BilinearFormEJA(B)
2153 sage: J1 = JordanSpinEJA(3)
2154 sage: J0.multiplication_table() == J0.multiplication_table()
2157 An error is raised if the matrix `B` does not correspond to a
2158 positive-definite bilinear form::
2160 sage: B = matrix.random(QQ,2,3)
2161 sage: J = BilinearFormEJA(B)
2162 Traceback (most recent call last):
2164 ValueError: bilinear form is not positive-definite
2165 sage: B = matrix.zero(QQ,3)
2166 sage: J = BilinearFormEJA(B)
2167 Traceback (most recent call last):
2169 ValueError: bilinear form is not positive-definite
2173 We can create a zero-dimensional algebra::
2175 sage: B = matrix.identity(AA,0)
2176 sage: J = BilinearFormEJA(B)
2180 We can check the multiplication condition given in the Jordan, von
2181 Neumann, and Wigner paper (and also discussed on my "On the
2182 symmetry..." paper). Note that this relies heavily on the standard
2183 choice of basis, as does anything utilizing the bilinear form matrix::
2185 sage: set_random_seed()
2186 sage: n = ZZ.random_element(5)
2187 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2188 sage: B11 = matrix.identity(QQ,1)
2189 sage: B22 = M.transpose()*M
2190 sage: B = block_matrix(2,2,[ [B11,0 ],
2192 sage: J = BilinearFormEJA(B)
2193 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2194 sage: V = J.vector_space()
2195 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2196 ....: for ei in eis ]
2197 sage: actual = [ sis[i]*sis[j]
2198 ....: for i in range(n-1)
2199 ....: for j in range(n-1) ]
2200 sage: expected = [ J.one() if i == j else J.zero()
2201 ....: for i in range(n-1)
2202 ....: for j in range(n-1) ]
2203 sage: actual == expected
2206 def __init__(self
, B
, field
=AA
, **kwargs
):
2209 if not B
.is_positive_definite():
2210 raise ValueError("bilinear form is not positive-definite")
2212 V
= VectorSpace(field
, n
)
2213 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2222 z0
= (B
*x
).inner_product(y
)
2223 zbar
= y0
*xbar
+ x0
*ybar
2224 z
= V([z0
] + zbar
.list())
2225 mult_table
[i
][j
] = z
2227 # Inner products are real numbers and not algebra
2228 # elements, so once we turn the algebra element
2229 # into a vector in inner_product(), we never go
2230 # back. As a result -- contrary to what we do with
2231 # self._multiplication_table -- we store the inner
2232 # product table as a plain old matrix and not as
2233 # an algebra operator.
2235 self
._inner
_product
_matrix
= ip_table
2237 super(BilinearFormEJA
, self
).__init
__(field
,
2242 # The rank of this algebra is two, unless we're in a
2243 # one-dimensional ambient space (because the rank is bounded
2244 # by the ambient dimension).
2245 self
.rank
.set_cache(min(n
,2))
2248 self
.one
.set_cache( self
.zero() )
2250 self
.one
.set_cache( self
.monomial(0) )
2253 def _max_random_instance_size():
2255 The maximum dimension of a random BilinearFormEJA.
2260 def random_instance(cls
, field
=AA
, **kwargs
):
2262 Return a random instance of this algebra.
2264 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2266 B
= matrix
.identity(field
, n
)
2267 return cls(B
, field
, **kwargs
)
2269 B11
= matrix
.identity(field
,1)
2270 M
= matrix
.random(field
, n
-1)
2271 I
= matrix
.identity(field
, n
-1)
2272 alpha
= field
.zero()
2273 while alpha
.is_zero():
2274 alpha
= field
.random_element().abs()
2275 B22
= M
.transpose()*M
+ alpha
*I
2277 from sage
.matrix
.special
import block_matrix
2278 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2281 return cls(B
, field
, **kwargs
)
2284 class JordanSpinEJA(BilinearFormEJA
):
2286 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2287 with the usual inner product and jordan product ``x*y =
2288 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2293 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2297 This multiplication table can be verified by hand::
2299 sage: J = JordanSpinEJA(4)
2300 sage: e0,e1,e2,e3 = J.gens()
2316 We can change the generator prefix::
2318 sage: JordanSpinEJA(2, prefix='B').gens()
2323 Ensure that we have the usual inner product on `R^n`::
2325 sage: set_random_seed()
2326 sage: J = JordanSpinEJA.random_instance()
2327 sage: x,y = J.random_elements(2)
2328 sage: actual = x.inner_product(y)
2329 sage: expected = x.to_vector().inner_product(y.to_vector())
2330 sage: actual == expected
2334 def __init__(self
, n
, field
=AA
, **kwargs
):
2335 # This is a special case of the BilinearFormEJA with the identity
2336 # matrix as its bilinear form.
2337 B
= matrix
.identity(field
, n
)
2338 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2341 def _max_random_instance_size():
2343 The maximum dimension of a random JordanSpinEJA.
2348 def random_instance(cls
, field
=AA
, **kwargs
):
2350 Return a random instance of this type of algebra.
2352 Needed here to override the implementation for ``BilinearFormEJA``.
2354 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2355 return cls(n
, field
, **kwargs
)
2358 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
,
2359 ConcreteEuclideanJordanAlgebra
):
2361 The trivial Euclidean Jordan algebra consisting of only a zero element.
2365 sage: from mjo.eja.eja_algebra import TrivialEJA
2369 sage: J = TrivialEJA()
2376 sage: 7*J.one()*12*J.one()
2378 sage: J.one().inner_product(J.one())
2380 sage: J.one().norm()
2382 sage: J.one().subalgebra_generated_by()
2383 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2388 def __init__(self
, field
=AA
, **kwargs
):
2390 self
._inner
_product
_matrix
= matrix(field
,0)
2391 super(TrivialEJA
, self
).__init
__(field
,
2395 # The rank is zero using my definition, namely the dimension of the
2396 # largest subalgebra generated by any element.
2397 self
.rank
.set_cache(0)
2398 self
.one
.set_cache( self
.zero() )
2401 def random_instance(cls
, field
=AA
, **kwargs
):
2402 # We don't take a "size" argument so the superclass method is
2403 # inappropriate for us.
2404 return cls(field
, **kwargs
)
2406 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2408 The external (orthogonal) direct sum of two other Euclidean Jordan
2409 algebras. Essentially the Cartesian product of its two factors.
2410 Every Euclidean Jordan algebra decomposes into an orthogonal
2411 direct sum of simple Euclidean Jordan algebras, so no generality
2412 is lost by providing only this construction.
2416 sage: from mjo.eja.eja_algebra import (random_eja,
2418 ....: RealSymmetricEJA,
2423 sage: J1 = HadamardEJA(2)
2424 sage: J2 = RealSymmetricEJA(3)
2425 sage: J = DirectSumEJA(J1,J2)
2433 The external direct sum construction is only valid when the two factors
2434 have the same base ring; an error is raised otherwise::
2436 sage: set_random_seed()
2437 sage: J1 = random_eja(AA)
2438 sage: J2 = random_eja(QQ)
2439 sage: J = DirectSumEJA(J1,J2)
2440 Traceback (most recent call last):
2442 ValueError: algebras must share the same base field
2445 def __init__(self
, J1
, J2
, **kwargs
):
2446 if J1
.base_ring() != J2
.base_ring():
2447 raise ValueError("algebras must share the same base field")
2448 field
= J1
.base_ring()
2450 self
._factors
= (J1
, J2
)
2454 V
= VectorSpace(field
, n
)
2455 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2459 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2460 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2464 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2465 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2467 super(DirectSumEJA
, self
).__init
__(field
,
2471 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2476 Return the pair of this algebra's factors.
2480 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2481 ....: JordanSpinEJA,
2486 sage: J1 = HadamardEJA(2,QQ)
2487 sage: J2 = JordanSpinEJA(3,QQ)
2488 sage: J = DirectSumEJA(J1,J2)
2490 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2491 Euclidean Jordan algebra of dimension 3 over Rational Field)
2494 return self
._factors
2496 def projections(self
):
2498 Return a pair of projections onto this algebra's factors.
2502 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2503 ....: ComplexHermitianEJA,
2508 sage: J1 = JordanSpinEJA(2)
2509 sage: J2 = ComplexHermitianEJA(2)
2510 sage: J = DirectSumEJA(J1,J2)
2511 sage: (pi_left, pi_right) = J.projections()
2512 sage: J.one().to_vector()
2514 sage: pi_left(J.one()).to_vector()
2516 sage: pi_right(J.one()).to_vector()
2520 (J1
,J2
) = self
.factors()
2523 V_basis
= self
.vector_space().basis()
2524 # Need to specify the dimensions explicitly so that we don't
2525 # wind up with a zero-by-zero matrix when we want e.g. a
2526 # zero-by-two matrix (important for composing things).
2527 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2528 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2529 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2530 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2531 return (pi_left
, pi_right
)
2533 def inclusions(self
):
2535 Return the pair of inclusion maps from our factors into us.
2539 sage: from mjo.eja.eja_algebra import (random_eja,
2540 ....: JordanSpinEJA,
2541 ....: RealSymmetricEJA,
2546 sage: J1 = JordanSpinEJA(3)
2547 sage: J2 = RealSymmetricEJA(2)
2548 sage: J = DirectSumEJA(J1,J2)
2549 sage: (iota_left, iota_right) = J.inclusions()
2550 sage: iota_left(J1.zero()) == J.zero()
2552 sage: iota_right(J2.zero()) == J.zero()
2554 sage: J1.one().to_vector()
2556 sage: iota_left(J1.one()).to_vector()
2558 sage: J2.one().to_vector()
2560 sage: iota_right(J2.one()).to_vector()
2562 sage: J.one().to_vector()
2567 Composing a projection with the corresponding inclusion should
2568 produce the identity map, and mismatching them should produce
2571 sage: set_random_seed()
2572 sage: J1 = random_eja()
2573 sage: J2 = random_eja()
2574 sage: J = DirectSumEJA(J1,J2)
2575 sage: (iota_left, iota_right) = J.inclusions()
2576 sage: (pi_left, pi_right) = J.projections()
2577 sage: pi_left*iota_left == J1.one().operator()
2579 sage: pi_right*iota_right == J2.one().operator()
2581 sage: (pi_left*iota_right).is_zero()
2583 sage: (pi_right*iota_left).is_zero()
2587 (J1
,J2
) = self
.factors()
2590 V_basis
= self
.vector_space().basis()
2591 # Need to specify the dimensions explicitly so that we don't
2592 # wind up with a zero-by-zero matrix when we want e.g. a
2593 # two-by-zero matrix (important for composing things).
2594 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2595 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2596 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2597 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2598 return (iota_left
, iota_right
)
2600 def inner_product(self
, x
, y
):
2602 The standard Cartesian inner-product.
2604 We project ``x`` and ``y`` onto our factors, and add up the
2605 inner-products from the subalgebras.
2610 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2611 ....: QuaternionHermitianEJA,
2616 sage: J1 = HadamardEJA(3,QQ)
2617 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2618 sage: J = DirectSumEJA(J1,J2)
2623 sage: x1.inner_product(x2)
2625 sage: y1.inner_product(y2)
2627 sage: J.one().inner_product(J.one())
2631 (pi_left
, pi_right
) = self
.projections()
2637 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2641 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance